0
# Integrationsgewichte (später hinterfragen...)
me.G5 = read.csv('xy.csv',sep=';',header=FALSE)[,2:3] 

me.x = me.G5[,1] 
me.w = me.G5[,2]
int.x = c((me.x/2+.5)*0.1,0.1+(me.x/2+.5)*0.9,1+(me.x/2+.5)*9,10+(me.x/2+.5)*90,100+(me.x/2+.5)*900)
int.w = c(me.w*0.1/2,me.w*.9/2,me.w*9/2,me.w*90/2,me.w*900/2)


### pdf
me.pdf  = function(MF,l) {
    pdf = function(x){
            for(i in 1:length(l)) X[,i] = MF[[i]](x)
            exp(l%*%t(X))
        }
    pdf
}

此代码在 R .... 我知道我们只使用列 2:3 。所以我理解代码,这是我的情节。

在此处输入图像描述

现在我用 Python 做了,我知道有很多库可以用于 pdf,我做到了

data =    np.loadtxt('xy.csv',skiprows=1, delimiter=';')

x =data[:,1]
w = data[:,2]

#int.x
a = (x/2 + 0.5)* 0.1 
b= 0.1+(x/2+.5)*0.9
g = 1+(x/2+.5)*9
c= 10+(x/2+.5)*90
d= 100+(x/2+.5)*900

X=[]
X.append(a)
X.append(b)
X.append(b2)
X.append(c)
X.append(d)

#int.w
 q= w* 0.1/2
e= w*0.9/2
r= w*9/2
m= w*90/2
n = w*900/2

W=[]
W.append(q)
W.append(e)
W.append(r)
W.append(m)
W.append(n)


plt.plot(X, W, marker='.', linestyle='none')
plt.margins(0.02)

我的做法很糟糕,我知道。我真的提高了我在 python 方面的技能。

............

在此处输入图像描述

这部分只是 R 代码的前 5 行:/。所以现在,我使用了 pdf 的库并让我成为了一个 norm-pdf 来比较两者。………………

from scipy.stats import norm
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1)
a= np.linspace(norm.ppf(0.01),norm.ppf(0.99),100)
ax.plot(a, norm.pdf(a), 'r-', lw=5, alpha=0.6, label='norm pdf')

这是这段代码的图表(norm-pdf) 在此处输入图像描述

所以我尝试做的就像R中的第一张图一样。或者它不能用于R中的定义pdf是Python中的库吗?!

4

1 回答 1

0

好吧对不起,我忘记了

这是我在 R 中的全部代码

### l
me.l    = function(MF,mt,sv=NA,tol=1e-10) {
    V   = mt[2:length(mt)]
    if(is.na(sv[1])) l  = rep(0,length(V)) else l = sv[-1]      
    X   = matrix(NaN,length(int.x),length(l))
    for(i in 2:length(MF)) X[,i-1]  = MF[[i]](int.x)
    run = 0;sig = 1
    while((max(abs(sig))>tol)&(run<40)) {
            m1  = as.vector(exp(l%*%t(X))*int.w)
            g   = m1%*%t(t(X)-V)
            G   = (t(X)-V)%*%(m1*t(t(X)-V))
            sig = solve(G,-t(g))
            l   = l + sig
            l   = as.vector(l)
            run = run +1            
        }
    m1  = as.vector(exp(l%*%t(X))*int.w)
    l0  = -log(sum(m1))
    c(l0,l) 
}



> MPO   = list();
> MPO[[1]] = function(x) 1
> MPO[[2]] = function(x) x
> MPO[[3]] = function(x) x^2
> MPO[[4]] = function(x) x^3
> l1 =  me.l(MPO[1:3],c(1,5,26))
> x = seq(0,10,,1e4)
> plot(x,me.pdf(MPO[1:3],l1)(x),type='l',lwd=2)
> lines(x,dnorm(x,5,1),col='red')

最后一部分是我在 R 控制台上写的。

现在我做到了

from scipy.stats import norm
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1)
a= np.linspace(norm.ppf(0.01),norm.ppf(0.99),100)
ax.plot(a, norm.pdf(a), 'r-', lw=5, alpha=0.6, label='norm pdf')
ax.plot(W, norm.pdf(W), '.')

我得到这张图 在此处输入图像描述

于 2018-02-07T19:17:02.410 回答