*! version 1.1 21oct2009
*! computes mean,sd for censored normal creation
*! austinnichols@gmail.com
prog cnormp, rclass
version 10
syntax [anything] [, Censpoint(real 0) ]
tokenize `anything'
if real(`"`1'"')==. | `"`1'"'=="" error 198
if real(`"`2'"')<0 |real(`"`2'"')==. | `"`2'"'=="" error 198
tempname b
di as txt _n  "To produce a censored normal dist. with"
di "mean `1' and SD `2' and censoring"
di "point `censpoint', generate using an underlying"
di "normal distribution with"
mat `b'=(`1',`2')
qui mata: fobj("`b'")
return scalar s=`b'[1,2]
return scalar m=`b'[1,1]
di "mean " as res `b'[1,1] as txt " and SD " as res `b'[1,2]
end
version 10
mata:
void obj(todo,b,crit,g,H)
{
external x,y,tm,ts,b,c
m=b[1,1]
s=b[1,2]
F=normal((c-m)/s)
f=normalden((c-m)/s)
h=f/(1-F)
x=F*c+(1-F)*(m+h*s)
y=sqrt(F*(x-c)^2+(1-F)*s^2*(1-h*(h+(m-c)/s))+(1-F)*(m+h*s-x)^2)
crit=(tm-x)^2+(ts-y)^2
}
void fobj(string scalar beta)
{
external x,y,tm,ts,b,c
c=strtoreal(st_local("censpoint"))
init=st_matrix(beta)
tm=init[1,1]
ts=init[1,2]
S=optimize_init()
optimize_init_evaluator(S, &obj())
optimize_init_which(S,"min")
optimize_init_evaluatortype(S,"d0")
optimize_init_params(S,init)
p=optimize(S)
st_replacematrix(beta,p)
}
end
exit