# PUF model and probability distributions

We use the temperature-dependent model for SRAM-PUF as defined by Roel Maes. 
\begin{align*}
r_i^{(j)}=\left\{\begin{matrix}
0 , &\text{if } m_i+d_i\cdot t^{(j)} +w\leq \tau,\\
1 , &\text{if } m_i+d_i \cdot t^{(j)} +w> \tau.
\end{matrix}\right.
\end{align*}
With process variable $m_i \sim \mathcal{N}(\mu_M,\sigma_M^2)$, temperature dependence variable $d_i \sim \mathcal{N}(0,\sigma_D^2)$, noise variable $w_i^{(j)} \sim \mathcal{N}(0,\sigma_N^2)$, and temperature $t^{(j)}$ at time $j$.
Roel also defines some parameters $\lambda_1 = \sigma_N/\sigma_M$, $\lambda_2 = (\tau-\mu_M)/\sigma_M$, and $\theta= \sigma_N / \sigma_D $. 

We refer to the one-probability of cell $i$ as $\xi_i$ and in case of known state-variables, this is given as $\xi_i = \Pr(X=1) = 1- cdf_N(\tau - m_i-d_i \cdot t^{(j)} )=\Phi(\frac{-\tau + m_i + d_i \cdot t^{(j)} }{\sigma_N}) $.

It is our goal to estimate $\lambda_1,\lambda_2,\theta$ for our dataset.

We will use maximum-likelihood estimation to achieve this.

## Settings of parameters


```python
import numpy as np
from scipy.stats import norm 
import matplotlib.pyplot as plt
import srampufparest
# set the parameters
lambda12 = [0.1,0.004]
theta = 45;
```

# verify srampufparest.loadUnique()
Does it give the same histogram as the code that I used before ?


```python
Nsamples = 60
Ndevices = 96
datapath = r'C:\Users\lkusters\surfdrive\python\server_008104\20181125\Unique'
folder = datapath+'\\002_Temp_025'
_,merged1 = srampufparest.loadUnique(folder)
import os
data25 = [[] for i in range(Nsamples)]
for ii in range(Ndevices):
    count = 0;
    for fname in os.listdir(folder):
        if fname.startswith("Unique_dev%03d"%(ii+1) ):
            data =  np.fromfile(folder+'\\'+fname, np.uint8) 
            data_bits = np.unpackbits(data)
            data25[count].extend(data_bits)
            count += 1;
data25 = np.swapaxes(data25,0,1) # swap axes, s.t. 1st dim is cells, and second is observations
```

    Finished loading data for 96 devices, with 60 observations of 65536 cells
    returning observations [96,60,65536] and merged observations [6291456,60]



```python
# from srampufparest
from numpy.random import uniform
print(len(merged1))
N25 = [len(d) for d in merged1] # number of measurements per cell
sums25 = [sum(n) for n in merged1 ]
print(set(N25))
# now generate errorcount
errorcount25 = [sum25 if uniform(0,1,1)>(sum25/Nsamples) else Nsamples-sum25 for sum25 in sums25]
```

    6291456
    {60}



```python
# original method
from numpy.random import uniform
print(len(data25))
N25 = [len(d) for d in data25] # number of measurements per cell
sums25 = [sum(n) for n in data25 ]
print(set(N25))
# now generate errorcount
errorcount25 = [sum25 if uniform(0,1,1)>(sum25/Nsamples) else Nsamples-sum25 for sum25 in sums25]
```

    6291456
    {60}



```python
# ORIGINAL
# plot the histogram of the number of observed ones in each cell
from scipy.special import comb
pmfdata_sum25,edges = np.histogram(sums25,bins=Nsamples+1,range=(-.5,Nsamples+.5))
plt.plot(range(Nsamples+1),[count/len(data25) for count in pmfdata_sum25] ,'o', fillstyle='none')
plt.yscale('log')
plt.xlabel('k Number of ones for total %d measurements'%(Nsamples))
plt.ylabel('#cells with Number of ones = k ')
plt.title('$T=25^oC$')

# Functions srampuf
pmfhist,centers = srampufparest.getcounts1D(merged1)
plt.plot(centers,[count/len(merged1) for count in pmfhist] )

plt.show()
```


```python
print(centers)
```

    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]


# Verify fitting method

# Pdf of one-probability (no temperature)

We find that apparently the approximations by python (in integrator?) are not accurate enough, since the pdf does not sum to 1 (but instead 0.38). Therefore, we instead use the equations to approximate the cdf and then calculate the pdf by integrating (sum over small intervals). 

## Pdf of the one-probability
\begin{equation}
        p_0(\xi|\lambda_1,\lambda_2)
        = \frac{\lambda_1 \cdot \phi( \Phi^{-1}(\xi) \cdot \lambda_1  +\lambda_2) }{\phi(\Phi^{-1}(\xi))}
\end{equation}

First, we show that indeed the total probability does not sum to 1:

## Cdf of the one-probability


The PDF can be approximated as 
\begin{equation}
        p_0(\xi|\lambda_1,\lambda_2)
        = \Pr(\Xi \leq \xi + \Delta\xi) - \Pr(\Xi \leq \xi - \Delta\xi) ,
\end{equation}
with 
\begin{equation}
        \Pr(\Xi \leq \xi)
        = \Phi( \Phi^{-1}(\xi) \cdot \lambda_1  +\lambda_2) 
\end{equation}

Second show that we can use cdf and integrate to improve the approximation


```python
# First plot the pdf
l1 = lambda12[0]
l2 = lambda12[1]
NX = 1000 # number of steps

xx = [i for i in np.linspace(0,1,NX)]
yy= [norm.cdf(l1*norm.ppf(xi)+l2) for xi in xx]
pdf_one = np.diff(yy)
pdf_x = [x+0.5/NX for x in xx[:-1]]
plt.plot(pdf_x,pdf_one,'o', fillstyle='none')

NX = 1000 # number of steps

p0_p, p0_xi = srampufparest.pdfp0(l1,l2,NX) 

plt.plot(p0_xi,p0_p)
plt.yscale('log')
plt.xlabel('Pr(X=1)=xi')
plt.ylabel('Pr(Xi=xi)')
plt.show()
pred = sum(p0_p) # should be 1
print(pred, 'should be equal to ONE ')
```


```python
# Then verify that we can reproduce the histogram of the one counts
# plot the histogram of the number of observed ones in each cell
from scipy.special import comb
pmfhist,centers = srampufparest.getcounts1D(merged1)
plt.plot(centers,[count/len(merged1) for count in pmfhist] ,'o', fillstyle='none')
plt.yscale('log')
plt.xlabel('k Number of ones for total %d measurements'%(Nsamples))
plt.ylabel('#cells with Number of ones = k ')
plt.title('$T=25^oC$')

# Now the fit
l1 = lambda12[0]
l2 = lambda12[1]
NX = 10000 # number of steps

p0_p, p0_xi = srampufparest.pdfp0(l1,l2,NX) 

pmf_ones = [0]*(Nsamples+1)
for nones in range(Nsamples+1): # loop over number of ones
    pred = [(xi**nones)*((1-xi)**(Nsamples-nones))*pdfp for (xi,pdfp) in zip(p0_xi,p0_p) ]
    pmf_ones[nones] = comb(Nsamples,nones)*sum(pred)
        
plt.plot(range(Nsamples+1),pmf_ones)

plt.show()
pred = sum(p0_p) # should be 1
print(pred, 'should be equal to ONE ')
```


```python
# Now the temperature dependent behavior
# settings from Roel: 0.12,0.02,45
l1 = 0.12
l2 = 0.02
theta = 45
NX = 1000 # number of steps

idx = srampufparest.getcellskones(merged2,5); # select cells that observed 5 ones in T2
merged1_sel = [merged1[i] for i in idx] # find observations of these cells at T1
hist,centerpoints = srampufparest.getcounts1D(merged1_sel); # generate the histogram

p1 = 5/60 # one probability
L = 40 # number of observations at T1
dT = 65 # temperature difference between T1 and T2
predict,x = srampufparest.pmfpkones_temperature(NX,dT,l1,l2,theta,p1,L); # predict the histogram

plt.plot(x,predict)
plt.plot(centerpoints,hist,'o', fillstyle='none') # actual values
plt.yscale('log')
plt.xlabel('k Number of ones for total %d measurements'%(Nsamples))
plt.ylabel('#cells with Number of ones = k ')

```

    [0.31010513237186754, 0.10058079065038165, 0.021414246053697624, 0.007509174365105543, 0.005113879303033711, 0.004256694020444828, 0.003723684383966019, 0.0033391881623195127, 0.0030468769351748547, 0.002817159259993022, 0.0026321161129814253, 0.00248018010332088, 0.0023535354726852997, 0.00224670598903209, 0.002155738815462368, 0.0020777093525364064, 0.0020104084943440115, 0.0019521382792435193, 0.001901574438489749, 0.0018576716049877391, 0.0018195965162510464, 0.0017866800605659812, 0.0017583823008481784, 0.0017342666268835474, 0.0017139804573404471, 0.0016972407335269588, 0.0016838229890183372, 0.001673553145528366, 0.0016663014385726656, 0.00166197805605192, 0.0016605302043487897, 0.0016619404168654263, 0.0016662260009514534, 0.0016734395896692502, 0.0016836708320020091, 0.0016970493256213913, 0.0017137489773863103, 0.0017339940771341318, 0.0017580675018859808, 0.0017863216473234089, 0.001819192936825101, 0.001857221125029318, 0.0019010751557642824, 0.0019515881560259484, 0.002009805420579908, 0.0020770512619482056, 0.002155023893422279, 0.0022459330418780135, 0.0023527045873945815, 0.002479293844995779, 0.002631181756850738, 0.002816193112185558, 0.0030459131712021817, 0.003338299040791898, 0.0037230470696340955, 0.004257062598871042, 0.005121381785408889, 0.007573918461233495, 0.021774823635995335, 0.10113958971143032, 0.30592430516568947]


# Pdf of one-probability (with temperature)


## Pdf of the one-probability
\begin{equation}
        p_1(\xi_2|\xi_1,t_1,t_2,\lambda_1,\lambda_2,\theta)
        = \phi \left(\frac{\theta}{t_2-t_1} \cdot (\Phi^{-1}(\xi_{t_2})-\Phi^{-1}(\xi_{t_1})) \right)   \cdot \frac{\theta}{t_2-t_1}  \cdot  \frac{1 }{\phi(\Phi^{-1}(\xi_{t_2} )  ) }
\end{equation}

First, we show that indeed the total probability does not sum to 1:

## Cdf of the one-probability


The PDF can be approximated as 
\begin{equation}
        p_1(\xi_2|\xi_1,t_1,t_2,\lambda_1,\lambda_2,\theta)
        = \Pr(\Xi_{t_2} \leq \xi_{t_2}+ \Delta\xi|\Xi_{t_1}=\xi_{t_1})  - \Pr(\Xi_{t_2} \leq \xi_{t_2}- \Delta\xi|\Xi_{t_1}=\xi_{t_1}) ,
\end{equation}
with 
\begin{equation}
        \Pr(\Xi_{t_2} \leq \xi_{t_2}|\Xi_{t_1}=\xi_{t_1}) 
        = \Phi \left(\frac{\theta}{t_2-t_1} \cdot (\Phi^{-1}(\xi_{t_2})-\Phi^{-1}(\xi_{t_1})) \right)
\end{equation}

We show that we can use cdf and integrate to improve the approximation


```python
# use the CDF
l1 = lambda12[0]
l2 = lambda12[1]
th = theta
NX = 100 # number of steps

dT = 25 # temperature difference

# first we choose single value for pi
pi = 0.1 # xit1
xx = [i for i in np.linspace(0,1,NX)]
yy = [norm.cdf( (th/dT) * (norm.ppf(xi)-norm.ppf(pi)) )   for xi in xx]
p1_p = np.diff(yy)
p1_xi = [x+0.5/NX for x in xx[:-1]]
plt.plot(p1_xi,p1_p)
plt.yscale('log')
plt.xlabel('Pr(X=1)=xi')
plt.ylabel('Pr(Xi=xi|pi=%0.2f,dT=%d)'%(pi,dT))
plt.show()
pred = sum(p1_p)
xx[0]

print(pred, 'should be equal to ONE ')

# second we integrate over all possible pi, to get the average distribution
# we re-use p0_p and p0_xi
ytotal = [ sum( [ Pdfpi*norm.cdf( (th/dT) * (norm.ppf(xi)-norm.ppf(pi)) )  for (Pdfpi,pi) in zip(p0_p,p0_xi) ])   for xi in xx]
ytotal = np.diff(ytotal)
ytotalx = [x+0.5/NX for x in xx[:-1]]

plt.plot(ytotalx,ytotal)
plt.yscale('log')
plt.xlabel('Pr(X=1)=xi')
plt.ylabel('Pr(Xi=xi|dT=%d)'%dT)
plt.show()
pred = sum(ytotal)

print(pred, 'should be equal to ONE ')
```


```python
# use the CDF
l1 = lambda12[0]
l2 = lambda12[1]
th = theta
NX = 100 # number of steps

dT = 25 # temperature difference

# first we choose single value for pi
pi = 0.1 # xit1
p1_p, p1_xi = srampufparest.pdfp1(pi,dT,l1,l2,th,NX)
plt.plot(p1_xi,p1_p)
plt.yscale('log')
plt.xlabel('Pr(X=1)=xi')
plt.ylabel('Pr(Xi=xi|pi=%0.2f,dT=%d)'%(pi,dT))
plt.show()
pred = sum(p1_p)
xx[0]

print(pred, 'should be equal to ONE ')

```


```python
# plot the histogram of the number of observed ones in each cell
pmfdata_sum25,edges = np.histogram(sums25,bins=Nsamples+1,range=(-.5,Nsamples+.5))
plt.plot(range(Nsamples+1),[count/len(data25) for count in pmfdata_sum25] ,'o', fillstyle='none')
plt.yscale('log')
plt.xlabel('k Number of ones for total %d measurements'%(Nsamples))
plt.ylabel('#cells with Number of ones = k ')
plt.title('$T=25^oC$')

# now also plot the fit
pmf_est = calculatePmf_kones_noT(lambda12,Nsamples)
plt.plot(range(Nsamples+1),pmf_est)
fig1 = plt.gcf()
plt.show()
fig1.savefig('UNIQUE_kones_pdf.jpg', bbox_inches='tight')
```

# Build the grid for ML estimation

Now the probability that for any cell we observe $k_1$ ones at temperature $t_1$ and $k_2$ ones at temperature $t_2$, is 
\begin{equation}
    \Pr_{k_1,k_2}(k_1,k_2 |t_1,t_2,\lambda_1,\lambda_2,\theta)
     ={K_1 \choose k_1} {K_2 \choose k_2}  \int_0^1  \xi_1^{k_1} (1-\xi_1)^{K_1 -k_1} p_0(\xi_1|\lambda_1,\lambda_2)   \int_0^1 \xi_2^{k_2} (1-\xi_2)^{K_2 -k_2} p_1(\xi_2|\xi_1,t_1,t_2 ,\lambda_1,\lambda_2,\theta) d \xi_2  d \xi_1 .
\end{equation}


```python

```


```python
from scipy.special import comb
# preparations
NX = 100
K=60 # number of observations
l1 = 0.1213;
l2 = 0.0210;
th = 45;
dT = 25 # difference between temperatures

xx = [i for i in np.linspace(0,1,NX)]
yy= [norm.cdf(l1*norm.ppf(xi)+l2) for xi in xx]
p0_p = np.diff(yy)
p0_xi = [x+0.5/NX for x in xx[:-1]]
```


```python
from scipy.special import comb
import sys
K=60

for th in [45,40,50,35,55]:
    print('Likelihood Pr_{k_1,k_2}(k_1,k_2 |dT,l1,l2,l2), with l1=%0.5f,l2=%0.5f,dT=%d,theta=%d'%(l1, l2, dT, th))

    Pk1k2 = [[0]*(K+1) for k in range(K+1)]
    for k1 in range(K+1):
        comb1 = comb(K,k1)
        for k2 in range(K+1):
            sys.stdout.write("\rCalculating k1=%2.0d,k2=%2.0d" % (k1,k2) )
            total = 0
            for (xi1,pdfxi1) in zip(p0_xi,p0_p):
                yy = [norm.cdf( (th/dT) * (norm.ppf(xi)-norm.ppf(xi1)) )   for xi in xx]
                p1_p = np.diff(yy)
                p1_xi = [x+0.5/NX for x in xx[:-1]]
                int2 = sum([(xi2**k2)*((1-xi2)**(K-k2))*pdfxi2 for (xi2,pdfxi2) in zip(p1_xi,p1_p)])
                total = total+ pdfxi1*int2*(xi1**k1)*((1-xi1)**(K-k1))
            Pk1k2[k1][k2] = comb1*comb(K,k2)*total
    np.savetxt('Pk1k1_(%d).txt.gz' %th , X, delimiter=' ', newline='\n', header='Likelihood Pr_{k_1,k_2}(k_1,k_2 |dT,l1,l2,l2), with l1=%0.5f,l2=%0.5f,dT=%d,theta=%d'%(l1, l2, dT, th), comments='# ')
```

    Likelihood Pr_{k_1,k_2}(k_1,k_2 |dT,l1,l2,l2), with l1=0.12130,l2=0.02100,dT=25,theta=45
    Calculating k1= 1,k2=50


```python
np.savetxt('Pk1k1_(%d).txt.gz' %th , X, delimiter=' ', newline='\n', header='Likelihood Pr_{k_1,k_2}(k_1,k_2 |dT,l1,l2,l2), with l1=%0.5f,l2=%0.5f,dT=%d,theta=%d'%(l1, l2, dT, th), comments='# ')
```

# Appendix : We can do the same for the pdf of Pe
Again we show that the pdf (approximate by python) does not sum to 1


```python
# can we also plot the pdf of Pe ?
lambda12 = [0.13,0.02]
l1 = lambda12[0]
l2 = lambda12[1]
NX = 1000 # number of steps

xx = [i for i in np.linspace(0,1,NX)]
y = [l1*(1-xi)*(norm.pdf(l1*norm.ppf(xi)+l2)+norm.pdf(l1*norm.ppf(xi)-l2))/norm.pdf(norm.ppf(xi)) for xi in xx]
plt.plot(xx,y)
plt.yscale('log')
plt.show()
pred = sum(y[1:-1])/NX
print(pred) # is not even close to 1 !
# this shows us that we are not accurate enough!
```

Again show that we can use cdf and integrate to improve the approximation


```python
# what if we use cdf instead?
lambda12 = [0.12,0.02]
l1 = lambda12[0]
l2 = lambda12[1]
NX = 100 # number of steps

xx = [i for i in np.linspace(0,1,NX)]
upperlim = [norm.ppf(i) for i in xx] # these are the upper limits of the integral I want to calculate
plt.plot(xx,upperlim)
plt.yscale('log');plt.xlabel('x');plt.ylabel(r'$\Phi^{-1}(x)$')
plt.title('Upper limit of integral')
plt.show()

# Now we generate the part inside the integral
Nu = 500000;
uu = np.linspace(-295,45,Nu) # outside of this range it is 0 anyway
y = [norm.cdf(-u)*(norm.pdf(l1*u+l2)+norm.pdf(l1*u-l2)) for u in uu];
plt.plot(uu,y)
#print(y[90:-355])
plt.yscale('log');
plt.xlabel('u');plt.ylabel(r'$\Phi(-u)\cdot(\phi(\dots)+\dots)$')
plt.title('Values of function inside integral')
plt.show()

# now I wil calculate the values
idxs = [np.searchsorted(uu, upper, side="left") for upper in upperlim];
cdfPe = [(340/Nu)*sum(y[:idx])*l1 for idx in idxs]
plt.plot(xx,cdfPe)
plt.yscale('log');plt.xlabel('x');plt.ylabel(r'$\lambda_1 \int_{-\infty}^{\Phi^{-1}(x)}\dots du$')
plt.title('$cdf_{Pe}$')
plt.show()

# ok! now I have cdf. Can I make pdf ?
pdfPe = np.diff(cdfPe)
plt.plot([x+0.5/NX for x in xx[:-1]],pdfPe)
plt.yscale('log');plt.xlabel('x');plt.ylabel('$\Delta cdf_{Pe}$')
plt.title('$pdf_{Pe}$')
plt.show()

print(sum(pdfPe))
```


```python

```

d512e8ef5ecc7f6e90f2a2aa06a407edd550c5d9
ipynb
Jupyter Notebook
srampufparest/tests/.ipynb_checkpoints/Verify module srampufparest-checkpoint.ipynb
lkusters/srampufparest
3c44b3c719db7071e17a5e830e5eecade16ce736
MIT
srampufparest/tests/.ipynb_checkpoints/Verify module srampufparest-checkpoint.ipynb
lkusters/srampufparest
3c44b3c719db7071e17a5e830e5eecade16ce736
MIT
srampufparest/tests/.ipynb_checkpoints/Verify module srampufparest-checkpoint.ipynb
lkusters/srampufparest
3c44b3c719db7071e17a5e830e5eecade16ce736
MIT
Qwen/Qwen-72B
1. YES
2. YES