source("cgamma.R")
dpearson4 <- function (x, alpha, m, nu, lambda){
abs(cgamma(m+nu/2*1i)/gamma(m))^2/(alpha*beta(m-1/2, 1/2))*(1+((x-lambda)/alpha)^2)^(-m)*exp(-nu*atan((x-lambda)/alpha))
}
params <- function(alpha, m, nu, lambda){
r <- 2*(m-1);
ave <- lambda - alpha*nu/r;
sig2 <- alpha^2*(r^2+nu^2)/(r^2*(r-1))
skew <- -4*nu/(r-2)*sqrt((r-1)/(r^2+nu^2))
kurt <- 3*(r-1)*((r+6)*(r^2+nu^2)-8*r^2)/((r-2)*(r-3)*(r^2+nu^2)) - 3
c(ave, sig2, skew, kurt)
}
norm2 <- function(x){(t(x)%*%x)[1]}
findparams <- function(m, s2, sk, k){
f <- function(x){
norm2(params(x[1],x[2],x[3],x[4])-c(m, s2, sk, k))
}
res <- nlm(f, c(sqrt(2), 4, 2, 0));#c(1, 5.3, 0, 1));
res$estimate
}
t=seq(-5, 5, 0.01);
p <- findparams(0, 1, 0, 0.05)
print(p);
alpha <- p[1]; m <- p[2]; nu <- p[3]; lambda <- p[4];
#print(params(alpha, m, nu, lambda))
plotpearson<-function(alpha, m, nu, lambda){
plot(t, dnorm(t), type="l", col="red", lwd=3, xlim=c(-5, 5), ylim=c(0, 0.65));
lines(t, dpearson4(t, alpha, m, nu, lambda), type="l", col="blue", lwd=3)
}
#lambda<-0; nu<-6; alpha<-sqrt(nu); m<-(nu+1)/2;
#print(params(alpha, m, nu, lambda))
#plot(t, dpearson4(t, alpha, m, nu, lambda), type="l")