-
Notifications
You must be signed in to change notification settings - Fork 0
/
gendpd.R
81 lines (61 loc) · 2.24 KB
/
gendpd.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
library("numDeriv")
library("nleqslv")
newton <- function(f1,f2,tol = 1e-7 , x0 = rep(1,2), N = 1000){
h = 1e-7
i = 1;
x1 = x0
p = list()
while (i <= N) {
f1value=try(f1(x0),silent=T)
f2value=try(f2(x0),silent=T)
if(class(f1value)=="try-error"){return("Start with new solution")}
if(class(f2value)=="try-error"){return("Start with new solution")}
j1=try(grad(f1,x0))
j2=try(grad(f2,x0))
if(class(j1)=="try-error"){return("Start with new solution")}
if(class(j2)=="try-error"){return("Start with new solution")}
j=rbind(j1,j2)
f=rbind(f1value,f2value)
if(det(j)==0){return("Start with new solution")}
invj=try(solve(j),silent=TRUE)
if(class(invj)=="try-error"){return("Start with new solution")}
x1 = (x0 - c(invj%*%f))
p[[i]] = x1
i = i + 1
if (norm(abs(x1 - x0),type="2") < tol) break
x0 = x1
}
return(p)
}
gendpd=function(smp,alp,beta,x0){
v=seq(0,2000,1)
fmd=function(p){if(p[2]>=0){
sum((dnorm(smp,mean=p[1],sd=p[2])^alp)*((smp-p[1])/(p[2]^2))*exp(beta*dnorm(smp,mean=p[1],sd=p[2])))/length(smp)}
else{
retrun("Start with new solution")
}
}
fwint=function(p){
if(p[2]>=0){
try(integrate(function(x){(dnorm(x,mean=p[1],sd=p[2])^(1+alp))*((x-p[1])/(p[2]^2))*exp(beta*dnorm(x,mean=p[1],sd=p[2]))}
,-Inf,Inf,rel.tol=1e-7)$value,silent=T)}
else{
return("Start with new solution")
}
}
if(class(fwint)=="try-error"){return("Start with new solution")}
fmu=function(p){try(fmd(p)-fwint(p),silent=T)}
fsd=function(p){if(p[2]>=0){sum((dnorm(smp,mean=p[1],sd=p[2])^alp)*(((smp-p[1])^2/(p[2]^2)-1)/p[2])*exp(beta*dnorm(smp,mean=p[1],sd=p[2])))/length(smp)}
else{return("Start with new solution")}}
fqint=function(p){
if(p[2]>=0){
try(integrate(function(x){(dnorm(x,mean=p[1],sd=p[2])^(1+alp))*(((x-p[1])^2/(p[2]^2)-1)/p[2])*exp(beta*dnorm(x,mean=p[1],sd=p[2]))},-Inf,Inf,rel.tol=1e-7)$value,silent=T)}
else{
return("Start with new solution")}
}
if(class(fqint)=="try-error"){return("Start with new solution")}
fsigma=function(p){try(fsd(p)-fqint(p),silent=T)}
f=function(p){c(fmu(p),fsigma(p))}
k=multiroot(f,start=x0)
return(k)
}