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")