0

我想知道是否有人可以帮助我将以下代码从 MatLab 翻译成 Python。该方程用于确定截断正态分布的 99% 置信区间。

function sigma = var_truncNormal( a, b, mu, sigma, data )

x1 = (a-mu)/sigma * normpdf( a, mu, sigma );
x2 = (b-mu)/sigma * normpdf( b, mu, sigma );

cx = normcdf( b, mu, sigma) - normcdf( a, mu, sigma );

yhat = var( data(data>(mu-3000)&data<(mu+3000)) );
sigma2 = yhat/((1+(x1-x2)/cx - ((x1-x2)/cx)^2));
sigma = sqrt( sigma2 );

return;


function ci99 = GetCI99( data )

mu = median( data );
sigma = std( data );
fprintf( 1, 'initial sigma = %.1f\n', sigma );
sigma = var_truncNormal( mu-3000, mu+3000, mu, sigma, data );
fprintf( 1, 'updated sigma = %.1f\n', sigma );
sigma = var_truncNormal( mu-3000, mu+3000, mu, sigma, data );
fprintf( 1, 'updated sigma = %.1f\n', sigma );

ci99 = 2*mu-norminv( 0.01, mu, sigma );

figure( 'visible', 'off' );
hist( data, 5000:200:20000 );
axis( [5000 35000 0 550] );

hold;
[n2, xx] = ksdensity( data, 'npoints', 100 );
plot( xx, n2*length(data)*200, 'r' );

hdl = plot( xx, normpdf( xx, mu, sigma )*length(data)*200, 'k' );
set( hdl, 'linewidth', 2 );

line( [ci99 ci99], [0 550] );

print( '-dpdf', 'testFigure' );
close;

return;

我将不胜感激任何帮助。

4

2 回答 2

1

正确的公式是：

def normpdf_python(x, mu, sigma):
   return 1/(sigma*np.sqrt(2*np.pi))*np.exp(-1*(x-mu)**2/ (2*sigma**2) )

您忘记了 (2*sigma**2) 中的括号

顺便说一句：norm.pdf(x,loc,scale) 结果相同

于 2015-09-15T09:37:44.447 回答
0

我想你会很高兴从 matlab 转到 python/numpy。你只需要一行一行地浏览它。

例如函数的第一行：

x1 = (a-mu)/sigma * normpdf(a, mu, sigma);

python中的normpdf为：1/(sigma*sqrt(2*pi)) exp(-1 (a-mu)**2/2*sigma**2)。

所以我们可以在python中定义一个小函数：

a = numpy.array([3,4,5,2,2,2,1,3,3,3,3])
mu = 3.1
sigma = .7
def normpdf_python(x, mu, sigma):
   return 1/(sigma*sqrt(2*pi))*exp(-1*(x-mu)**2/2*sigma**2)

x1 = (a-mu)/sigma*normpdf_python(a,mu,sigma)

Numpy/Scipy 有一套很深的统计包，你应该经常检查预先存在的函数来满足你的需要。在这种情况下http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm pdf(x, loc=0, scale=1) 很接近，但是可能还不够，因为它被定义为：

norm.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
于 2013-11-12T14:31:28.930 回答