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
}

This Code is in R .... I know that we only use the columns 2:3 . So i understand the Code and this is my Plot .

enter image description here

Now i did it in Python, I know that there is many of libraries which i can use for pdf, which i did

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)

The way how i did it is very bad , i know it . And i really improve my skills in python.

...........

enter image description here

This parts was only the first 5 lines from the R code :/ . So now , i used the library for pdf's und made me a norm-pdf to compare the both. ......... ........

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

And this is the graph for this code (norm-pdf) enter image description here

So what i try to do is like the same as in the first graph from R. Or is it not possible to use for the definition pdf in R a library in Python?!

Moritz
  • 47
  • 6
  • Your W is not a vector but a list. When you append to the list, you get a list where each element is each np.array, not a large array. To do that you need to use `np.concatenate()` – Felipe Gerard Feb 07 '18 at 18:44
  • thanks for the information , i did it like here : W=np.concatenate((q,e,r,m,n)) But how can i compare it with norm-pdf – Moritz Feb 07 '18 at 18:59
  • Try using ax.plot again with the other data to plot one on top of the other. Also, note you didn't include your R plotting code, so it's hard to know how you did it or what you want. – Felipe Gerard Feb 07 '18 at 19:03

1 Answers1

0

okay Sorry , i forget it

This is my whole Code in 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')

The last part is what i write on the R console.

Now i did it

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), '.')

and i get this graph enter image description here

Moritz
  • 47
  • 6