Решаем систему линейных алгебраических уравнений с python-пакетом scipy.linalg (не путать с numpy.linalg)

Проверка на нормальность в Scipy

Нормальный закон распределения является простым и удобным для дальнейшего исследования. Чтобы проверить имеет ли тот или иной атрибут нормальное распределение, можно воспользоваться двумя критериями Python-библиотеки с модулем . Модуль   поддерживает большой диапазон статистических функций, полный перечень которых представлен в официальной документации.

В основе проверки на “нормальность” лежит проверка гипотез. Нулевая гипотеза – данные распределены нормально, альтернативная гипотеза – данные не имеют нормального распределения.

Проведем первый критерий Шапиро-Уилк [], возвращающий значение вычисленной статистики и p-значение. В качестве критического значения в большинстве случаев берется 0.05. При p-значении меньше 0.05 мы вынуждены отклонить нулевую гипотезу.

Проверим распределение атрибута Rings, количество колец:

import scipy

stat, p = scipy.stats.shapiro(data) # тест Шапиро-Уилк
print('Statistics=%.3f, p-value=%.3f' % (stat, p))

alpha = 0.05
if p > alpha:
    print('Принять гипотезу о нормальности')
else:
    print('Отклонить гипотезу о нормальности')

В результате мы получили низкое p-значение и, следовательно, отклоняем нулевую гипотезу:

Statistics=0.931, p-value=0.000
Отклонить гипотезу о нормальности

Второй тест по критерию согласия Пирсона [], который тоже возвращает соответствующее значение статистики и p-значение:

stat, p = scipy.stats.normaltest(data) # Критерий согласия Пирсона
print('Statistics=%.3f, p-value=%.3f' % (stat, p))

alpha = 0.05
if p > alpha:
    print('Принять гипотезу о нормальности')
else:
    print('Отклонить гипотезу о нормальности')

Этот критерий также отвергает нулевую гипотезу о нормальности распределения колец у моллюсков, так как p-значение меньше 0.05:

Statistics=242.159, p-value=0.000
Отклонить гипотезу о нормальности

File Input / Output package:

Scipy, I/O package, has a wide range of functions for work with different files format which are Matlab, Arff, Wave, Matrix Market, IDL, NetCDF, TXT, CSV and binary format.

Let us take one file format Python SciPy example as which are regularly used in MatLab:

 import numpy as np
 from scipy import io as sio
 array = np.ones((4, 4))
 sio.savemat('example.mat', {'ar': array}) 
 data = sio.loadmat(‘example.mat', struct_as_record=True)
 data

Output:

array(,
           ,
           ,
           ])

Code Explanation

  • Line 1 & 2: Import the essential SciPy library in Python with I/O package and Numpy.
  • Line 3: Create 4 x 4, dimensional one’s array
  • Line 4: Store array in example.mat file.
  • Line 5: Get data from example.mat file
  • Line 6: Print output.

Special Functions of Python SciPy

The following is the list of some of the most commonly used Special functions from the package of SciPy:

  • Cubic Root
  • Exponential Function
  • Log-Sum Exponential Function
  • Gamma

1. Cubic Root

The function is used to provide the element-wise cube root of the inputs provided.

Example:

from scipy.special import cbrt
val = cbrt()
print(val)

Output:

2. Exponential Function

The function is used to calculate the element-wise exponent of the given inputs.

Example:

from scipy.special import exp10
val = exp10()
print(val)

Output:

3. Log-Sum Exponential Function

The function is used to calculate the logarithmic value of the sum of the exponents of the input elements.

Example:

from scipy.special import logsumexp
import numpy as np
inp = np.arange(5)
val = logsumexp(inp)
print(val)

Here, numpy.arange() function is used to generate a sequence of numbers to be passed as input.

Output:

4.451914395937593

4. Gamma Function

Gamma function is used to calculate the gamma value, referred to as generalized factorial because, gamma(n+1) = n!

The function is used to calculate the gamma value of the input element.

Example:

from scipy.special import gamma
val = gamma()
print(val)

Output:

1.6.6. Statistics and random numbers: scipy.stats¶

The module contains statistical tools and probabilistic
descriptions of random processes. Random number generators for various
random process can be found in .

1.6.6.1. Distributions: histogram and probability density function

Given observations of a random process, their histogram is an estimator of
the random process’s PDF (probability density function):

>>> samples = np.random.normal(size=1000)
>>> bins = np.arange(-4, 5)
>>> bins
array()
>>> histogram = np.histogram(samples, bins=bins, density=True) + bins)
>>> bins
array()
>>> from scipy import stats
>>> pdf = stats.norm.pdf(bins)  # norm is a distribution object

>>> plt.plot(bins, histogram) 

>>> plt.plot(bins, pdf) 

The distribution objects

is a distribution object: each distribution
in is represented as an object. Here it’s the
normal distribution, and it comes with a PDF, a CDF, and much more.

If we know that the random process belongs to a given family of random
processes, such as normal processes, we can do a maximum-likelihood fit
of the observations to estimate the parameters of the underlying
distribution. Here we fit a normal process to the observed data:

>>> loc, std = stats.norm.fit(samples)
>>> loc     
-0.045256707...
>>> std     
0.9870331586...

Exercise: Probability distributions

Generate 1000 random variates from a gamma distribution with a shape
parameter of 1, then plot a histogram from those samples. Can you plot the
pdf on top (it should match)?

Extra: the distributions have many useful methods. Explore them by
reading the docstring or by using tab completion. Can you recover
the shape parameter 1 by using the method on your random
variates?

1.6.6.2. Mean, median and percentiles

The mean is an estimator of the center of the distribution:

>>> np.mean(samples)     
-0.0452567074...

The median another estimator of the center. It is the value with half of
the observations below, and half above:

>>> np.median(samples)     
-0.0580280347...

Tip

Unlike the mean, the median is not sensitive to the tails of the
distribution. It is “robust”.

Exercise: Compare mean and median on samples of a Gamma distribution

The median is also the percentile 50, because 50% of the observation are
below it:

>>> stats.scoreatpercentile(samples, 50)     
-0.0580280347...

Similarly, we can calculate the percentile 90:

>>> stats.scoreatpercentile(samples, 90)     
1.2315935511...

Tip

The percentile is an estimator of the CDF: cumulative distribution
function.

Metric (SI) Prefixes:

Return the specified unit in meter (e.g.
returns )

Example

from scipy import constants
print(constants.yotta)    #1e+24
print(constants.zetta)   
#1e+21print(constants.exa)     
#1e+18
print(constants.peta)     #1000000000000000.0print(constants.tera)    
#1000000000000.0
print(constants.giga)     #1000000000.0print(constants.mega)    
#1000000.0
print(constants.kilo)     #1000.0print(constants.hecto)   
#100.0
print(constants.deka)     #10.0print(constants.deci)    
#0.1
print(constants.centi)    #0.01print(constants.milli)   
#0.001
print(constants.micro)    #1e-06print(constants.nano)    
#1e-09
print(constants.pico)     #1e-12print(constants.femto)    #1e-15
print(constants.atto)     #1e-18print(constants.zepto)    #1e-21

SciPy (the library)¶

There is now a journal article available for citing usage of SciPy:

Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland,
Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson,
Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt,
Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov,
Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey,
İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde,
Josef Perktold, Robert Cimrman, Ian Henriksen, E.A. Quintero, Charles R Harris,
Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt,
and SciPy 1.0 Contributors. (2020) SciPy 1.0: Fundamental Algorithms
for Scientific Computing in Python. Nature Methods, 17(3), 261-272.

Here’s an example of a BibTeX entry:

@ARTICLE{2020SciPy-NMeth,
  author  = {Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E. and
            Haberland, Matt and Reddy, Tyler and Cournapeau, David and
            Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and
            Bright, Jonathan and {van der Walt}, St{\'e}fan J. and
            Brett, Matthew and Wilson, Joshua and Millman, K. Jarrod and
            Mayorov, Nikolay and Nelson, Andrew R. J. and Jones, Eric and
            Kern, Robert and Larson, Eric and Carey, C J and
            Polat, {\.I}lhan and Feng, Yu and Moore, Eric W. and
            {VanderPlas}, Jake and Laxalde, Denis and Perktold, Josef and
            Cimrman, Robert and Henriksen, Ian and Quintero, E. A. and
            Harris, Charles R. and Archibald, Anne M. and
            Ribeiro, Ant{\^o}nio H. and Pedregosa, Fabian and
            {van Mulbregt}, Paul and {SciPy 1.0 Contributors}},
  title   = {{{SciPy} 1.0 Fundamental Algorithms for Scientific
            Computing in Python}},
  journal = {Nature Methods},
  year    = {2020},
  volume  = {17},
  pages   = {261--272},
  adsurl  = {https//rdcu.beb08Wh},
  doi     = {10.1038s41592-019-0686-2},
}

Почему SciPy?

SciPy предоставляет высокоуровневые команды и классы для управления данными и визуализации данных, что значительно увеличивает мощность интерактивного сеанса Python.

Помимо математических алгоритмов в SciPy, программисту доступно все, от классов, веб-подпрограмм и баз данных до параллельного программирования, что упрощает и ускоряет разработку сложных и специализированных приложений.

Поскольку SciPy имеет открытый исходный код, разработчики по всему миру могут вносить свой вклад в разработку дополнительных модулей, что очень полезно для научных приложений, использующих SciPy.

Почему это важно учитывать?

Это связано с тем, что каждая итерация соответствуетвычислительные (а иногда и не вычислительные, а фактические физические) затраты,

Этот тип сценария возникает, когда оптимизация проводится не с простой математической оценкой, а с помощью сложного, трудоемкого моделирования или затратных и трудоемких экспериментов.

Идем сложнее — многовариантная функция

Хотя мы рассмотрели все существенные аспекты решения стандартной задачи оптимизации в предыдущих разделах, пример состоял из простой аналитической функции с одной переменной.

Но это не обязательно так!

Методы SciPy работают с любой функцией Python — не обязательно с замкнутой, одномерной математической функцией.

Покажем пример с многозначной функцией.

SciPy Ndimage

The SciPy provides the ndimage (n-dimensional image) package, that contains the number of general image processing and analysis functions. Some of the most common tasks in image processing are as follows:

  • Basic manipulations − Cropping, flipping, rotating, etc.
  • Image filtering − Denoising, sharpening, etc.
  • Image segmentation − Labeling pixels corresponding to different objects
  • Classification
  • Feature extraction
  • Registration

Here are some examples in which we will apply some of these image processing techniques on the images:

First, let us import an image that is already included in the SciPy package:

import scipy.misc
import matplotlib.pyplot as plt
face = scipy.misc.face()#returns an image of raccoon
#display image using matplotlib
plt.imshow(face)
plt.show()

Output:

Crop image

import scipy.misc
import matplotlib.pyplot as plt
face = scipy.misc.face()#returns an image of raccoon
lx,ly,channels= face.shape
# Cropping
crop_face = face[int(lx/4):int(-lx/4), int(ly/4):int(-ly/4)]
plt.imshow(crop_face)
plt.show()

Output:

Rotate Image

from scipy import misc,ndimage
import matplotlib.pyplot as plt
face = misc.face()
rotate_face = ndimage.rotate(face, 180)
plt.imshow(rotate_face)
plt.show()


Output:

Blurring or Smoothing  Images

Here we will blur the original images using the Gaussian filter and see how to control the level of smoothness using the sigma parameter.

from scipy import ndimage,misc
import matplotlib.pyplot as plt

face = scipy.misc.face(gray=True)
blurred_face = ndimage.gaussian_filter(face, sigma=3)
very_blurred = ndimage.gaussian_filter(face, sigma=5)

plt.figure(figsize=(9, 3))
plt.subplot(131)
plt.imshow(face, cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(132)
plt.imshow(very_blurred, cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(133)
plt.imshow(blurred_face, cmap=plt.cm.gray)
plt.axis('off')


plt.subplots_adjust(wspace=0, hspace=0., top=0.99, bottom=0.01,
                    left=0.01, right=0.99)

plt.show()

Output:

The first image is the original image followed by the blurred images with different sigma values.

Sharpening images

Here we will blur the image using the Gaussian method mentioned above and then sharpen the image by adding intensity to each pixel of the blurred image.

import scipy
from scipy import ndimage
import matplotlib.pyplot as plt

f = scipy.misc.face(gray=True).astype(float)
blurred_f = ndimage.gaussian_filter(f, 3)

filter_blurred_f = ndimage.gaussian_filter(blurred_f, 1)

alpha = 30
sharpened = blurred_f + alpha * (blurred_f - filter_blurred_f)

plt.figure(figsize=(12, 4))

plt.subplot(131)
plt.imshow(f, cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(132)
plt.imshow(blurred_f, cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(133)
plt.imshow(sharpened, cmap=plt.cm.gray)
plt.axis('off')

plt.tight_layout()
plt.show()

Output:

Edge detection

Edge detection includes a variety of mathematical methods that aim at identifying points in a digital image at which the image brightness changes sharply or, more formally, has discontinuities. The points at which image brightness changes sharply are typically organized into a set of curved line segments termed edges.

Here is an example where we make a square figure and then find its edges:

mport numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

im = np.zeros((256, 256))
im = 1

im = ndimage.rotate(im, 15, mode='constant')
im = ndimage.gaussian_filter(im, 8)

sx = ndimage.sobel(im, axis=0, mode='constant')
sy = ndimage.sobel(im, axis=1, mode='constant')
sob = np.hypot(sx, sy)

plt.figure(figsize=(9,5))
plt.subplot(141)
plt.imshow(im)
plt.axis('off')
plt.title('square', fontsize=20)
plt.subplot(142)
plt.imshow(sob)
plt.axis('off')
plt.title('Sobel filter', fontsize=20)
plt.show()

Output:

Length:

Return the specified unit in meters (e.g. returns )

Example

from scipy import constants

print(constants.inch)              #0.0254
print(constants.foot)              #0.30479999999999996
print(constants.yard)              #0.9143999999999999
print(constants.mile)              #1609.3439999999998
print(constants.mil)               #2.5399999999999997e-05
print(constants.pt)                #0.00035277777777777776
print(constants.point)             #0.00035277777777777776
print(constants.survey_foot)       #0.3048006096012192
print(constants.survey_mile)       #1609.3472186944373
print(constants.nautical_mile)     #1852.0
print(constants.fermi)             #1e-15
print(constants.angstrom)          #1e-10
print(constants.micron)            #1e-06
print(constants.au)                #149597870691.0
print(constants.astronomical_unit) #149597870691.0
print(constants.light_year)        #9460730472580800.0
print(constants.parsec)            #3.0856775813057292e+16

1.6.10. Image manipulation: scipy.ndimage¶

provides manipulation of n-dimensional arrays as
images.

1.6.10.1. Geometrical transformations on images

Changing orientation, resolution, ..

>>> from scipy import misc  # Load an image
>>> face = misc.face(gray=True)

>>> from scipy import ndimage # Shift, roate and zoom it
>>> shifted_face = ndimage.shift(face, (50, 50))
>>> shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest')
>>> rotated_face = ndimage.rotate(face, 30)
>>> cropped_face = face50-50, 50-50
>>> zoomed_face = ndimage.zoom(face, 2)
>>> zoomed_face.shape
(1536, 2048)

>>> plt.subplot(151)    
<matplotlib.axes._subplots.AxesSubplot object at 0x...>

>>> plt.imshow(shifted_face, cmap=plt.cm.gray)    
<matplotlib.image.AxesImage object at 0x...>

>>> plt.axis('off')
(-0.5, 1023.5, 767.5, -0.5)

>>> # etc.

1.6.10.2. Image filtering

Generate a noisy face:

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> face = face  # crop out square on right
>>> import numpy as np
>>> noisy_face = np.copy(face).astype(np.float)
>>> noisy_face += face.std() * 0.5 * np.random.standard_normal(face.shape)

Apply a variety of filters on it:

>>> blurred_face = ndimage.gaussian_filter(noisy_face, sigma=3)
>>> median_face = ndimage.median_filter(noisy_face, size=5)
>>> from scipy import signal
>>> wiener_face = signal.wiener(noisy_face, (5, 5))

Other filters in and
can be applied to images.

Exercise

1.6.10.3. Mathematical morphology

Tip

Mathematical morphology stems from set
theory. It characterizes and transforms geometrical structures. Binary
(black and white) images, in particular, can be transformed using this
theory: the sets to be transformed are the sets of neighboring
non-zero-valued pixels. The theory was also extended to gray-valued
images.

Mathematical-morphology operations use a structuring element
in order to modify geometrical structures.

Let us first generate a structuring element:

>>> el = ndimage.generate_binary_structure(2, 1)
>>> el 
array(,
       ,
       ])
>>> el.astype(np.int)
array(,
       ,
       ])
  • Erosion

    >>> a = np.zeros((7, 7), dtype=np.int)
    >>> a16, 25 = 1
    >>> a
    array(,
           ,
           ,
           ,
           ,
           ,
           ])
    >>> ndimage.binary_erosion(a).astype(a.dtype)
    array(,
           ,
           ,
           ,
           ,
           ,
           ])
    >>> #Erosion removes objects smaller than the structure
    >>> ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
    array(,
           ,
           ,
           ,
           ,
           ,
           ])
    
  • Dilation

    >>> a = np.zeros((5, 5))
    >>> a2, 2 = 1
    >>> a
    array(,
           ,
           ,
           ,
           ])
    >>> ndimage.binary_dilation(a).astype(a.dtype)
    array(,
           ,
           ,
           ,
           ])
    
  • Opening

    >>> a = np.zeros((5, 5), dtype=np.int)
    >>> a14, 14 = 1
    >>> a4, 4 = 1
    >>> a
    array(,
           ,
           ,
           ,
           ])
    >>> # Opening removes small objects
    >>> ndimage.binary_opening(a, structure=np.ones((3, 3))).astype(np.int)
    array(,
           ,
           ,
           ,
           ])
    >>> # Opening can also smooth corners
    >>> ndimage.binary_opening(a).astype(np.int)
    array(,
           ,
           ,
           ,
           ])
    
  • Closing:

Exercise

An opening operation removes small structures, while a closing operation
fills small holes. Such operations can therefore be used to “clean” an
image.

>>> a = np.zeros((50, 50))
>>> a10-10, 10-10 = 1
>>> a += 0.25 * np.random.standard_normal(a.shape)
>>> mask = a>=0.5
>>> opened_mask = ndimage.binary_opening(mask)
>>> closed_mask = ndimage.binary_closing(opened_mask)

Exercise

For gray-valued images, eroding (resp. dilating) amounts to replacing
a pixel by the minimal (resp. maximal) value among pixels covered by the
structuring element centered on the pixel of interest.

>>> a = np.zeros((7, 7), dtype=np.int)
>>> a16, 16 = 3
>>> a4, 4 = 2; a2, 3 = 1
>>> a
array(,
       ,
       ,
       ,
       ,
       ,
       ])
>>> ndimage.grey_erosion(a, size=(3, 3))
array(,
       ,
       ,
       ,
       ,
       ,
       ])

Discrete Fourier Transform – scipy.fftpack

  • DFT is a mathematical technique which is used in converting spatial data into frequency data.
  • FFT (Fast Fourier Transformation) is an algorithm for computing DFT
  • FFT is applied to a multidimensional array.
  • Frequency defines the number of signal or wavelength in particular time period.

Example: Take a wave and show using Matplotlib library. we take simple periodic function example of sin(20 × 2πt)

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np 

#Frequency in terms of Hertz
fre  = 5 
#Sample rate
fre_samp = 50
t = np.linspace(0, 2, 2 * fre_samp, endpoint = False )
a = np.sin(fre  * 2 * np.pi * t)
figure, axis = plt.subplots()
axis.plot(t, a)
axis.set_xlabel ('Time (s)')
axis.set_ylabel ('Signal amplitude')
plt.show()

Output:

You can see this. Frequency is 5 Hz and its signal repeats in 1/5 seconds – it’s call as a particular time period.

Now let us use this sinusoid wave with the help of DFT application.

from scipy import fftpack

A = fftpack.fft(a)
frequency = fftpack.fftfreq(len(a)) * fre_samp
figure, axis = plt.subplots()

axis.stem(frequency, np.abs(A))
axis.set_xlabel('Frequency in Hz')
axis.set_ylabel('Frequency Spectrum Magnitude')
axis.set_xlim(-fre_samp / 2, fre_samp/ 2)
axis.set_ylim(-5, 110)
plt.show()

Output:

  • You can clearly see that output is a one-dimensional array.
  • Input containing complex values are zero except two points.
  • In DFT example we visualize the magnitude of the signal.

NumPy¶

There is now a journal article available for citing usage of NumPy:

Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf
Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor,
Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan
Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime
Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant,
Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi,
Christoph Gohlke & Travis E. Oliphant.
Array programming with NumPy, Nature, 585, 357–362 (2020),
DOI:10.1038/s41586-020-2649-2 (publisher link)

Here’s an example of a BibTeX entry:

Image Processing with SciPy – scipy.ndimage

  • scipy.ndimage is a submodule of SciPy which is mostly used for performing an image related operation
  • ndimage means the “n” dimensional image.
  • SciPy Image Processing provides Geometrics transformation (rotate, crop, flip), image filtering (sharp and de nosing), display image, image segmentation, classification and features extraction.
  • MISC Package in SciPy contains prebuilt images which can be used to perform image manipulation task

Example: Let’s take a geometric transformation example of images

from scipy import misc
from matplotlib import pyplot as plt
import numpy as np
#get face image of panda from misc package
panda = misc.face()
#plot or show image of face
plt.imshow( panda )
plt.show()

Output:

Now we Flip-down current image:

#Flip Down using scipy misc.face image  
flip_down = np.flipud(misc.face())
plt.imshow(flip_down)
plt.show()

Output:

Example: Rotation of Image using Scipy,

from scipy import ndimage, misc
from matplotlib import pyplot as plt
panda = misc.face()
#rotatation function of scipy for image – image rotated 135 degree
panda_rotate = ndimage.rotate(panda, 135)
plt.imshow(panda_rotate)
plt.show()

Output:

Integration with Scipy – Numerical Integration

  • When we integrate any function where analytically integrate is not possible, we need to turn for numerical integration
  • SciPy provides functionality to integrate function with numerical integration.
  • scipy.integrate library has single integration, double, triple, multiple, Gaussian quadrate, Romberg, Trapezoidal and Simpson’s rules.

Example: Now take an example of Single Integration

Here a is the upper limit and b is the lower limit

from scipy import integrate
# take f(x) function as f
f = lambda x : x**2
#single integration with a = 0 & b = 1  
integration = integrate.quad(f, 0 , 1)
print(integration)

Output:

(0.33333333333333337, 3.700743415417189e-15)

Here function returns two values, in which the first value is integration and second value is estimated error in integral.

Example: Now take an SciPy example of double integration. We find the double integration of the following equation,

from scipy import integrate
import numpy as np
#import square root function from math lib
from math import sqrt
# set  fuction f(x)
f = lambda x, y : 64 *x*y
# lower limit of second integral
p = lambda x : 0
# upper limit of first integral
q = lambda y : sqrt(1 - 2*y**2)
# perform double integration
integration = integrate.dblquad(f , 0 , 2/4,  p, q)
print(integration)

Output:

(3.0, 9.657432734515774e-14)

You have seen that above output as same previous one.

1.6.4. Interpolation: scipy.interpolate¶

is useful for fitting a function from experimental
data and thus evaluating points where no measure exists. The module is based
on the FITPACK Fortran subroutines.

By imagining experimental data close to a sine function:

>>> measured_time = np.linspace(, 1, 10)
>>> noise = (np.random.random(10)*2 - 1) * 1e-1
>>> measures = np.sin(2 * np.pi * measured_time) + noise

can build a linear interpolation
function:

>>> from scipy.interpolate import interp1d
>>> linear_interp = interp1d(measured_time, measures)

Then the result can be evaluated at the time of interest:

>>> interpolation_time = np.linspace(, 1, 50)
>>> linear_results = linear_interp(interpolation_time)

A cubic interpolation can also be selected by providing the optional
keyword argument:

>>> cubic_interp = interp1d(measured_time, measures, kind='cubic')
>>> cubic_results = cubic_interp(interpolation_time)

1.6.9. Signal processing: scipy.signal¶

Tip

is for typical signal processing: 1D,
regularly-sampled signals.

Resampling : resample a signal to n
points using FFT.

>>> t = np.linspace(, 5, 100)
>>> x = np.sin(t)

>>> from scipy import signal
>>> x_resampled = signal.resample(x, 25)

>>> plt.plot(t, x) 

>>> plt.plot(t, x_resampled, 'ko') 

Tip

Notice how on the side of the window the resampling is less accurate
and has a rippling effect.

This resampling is different from the provided by as it
only applies to regularly sampled data.

Detrending : remove linear trend from signal:

>>> t = np.linspace(, 5, 100)
>>> x = t + np.random.normal(size=100)

>>> from scipy import signal
>>> x_detrended = signal.detrend(x)

>>> plt.plot(t, x) 

>>> plt.plot(t, x_detrended) 

Filtering:
For non-linear filtering, has filtering (median
filter , Wiener ),
but we will discuss this in the image section.

Tip

also has a full-blown set of tools for the design
of linear filter (finite and infinite response filters), but this is
out of the scope of this tutorial.

Spectral analysis:
compute a spectrogram –frequency
spectrums over consecutive time windows–, while
comptes a power spectrum density (PSD).

Sub-Packages in Python SciPy

There are various sub-modules available in the SciPy library to perform and enhance the efficiency of the scientific calculations.

Some of the popular sub-modules of the SciPy library are listed below:

  • special: This sub-module contains the Special functions to perform a specific task.
  • constants: Represents constants.
  • optimize: This sub-module contains algorithms for optimization.
  • integrate: This sub-module contains functions to perform Mathematical Integration.
  • interpolate: Represents functions to perform interpolation.
  • linalg: Represents functions to perform operations on linear algebra equations.
  • io: It contains functions to perform Input/Output operations on the given input.
  • fftpack: Represents functions to perform Discrete Fourier Transform.
  • signal: Represents functions and tools for Signal Processing in Python.
  • sparse: Represents algorithms to deal with sparse matrices.
  • cluster: Represents functions to perform hierarchical clustering.

SciPy ODR

Orthogonal Distance Regression (ODR) is the name given to the computational problem associated with finding the maximum likelihood estimators of parameters in measurement error models in the case of normally distributed errors.

The Least square method calculates the error vertical to the line (shown by grey colour here) whereas ODR calculates the error perpendicular(orthogonal) to the line. This accounts for the error in both X and Y whereas using  Least square method, we only consider the error in Y.

scipy.odr Implementation for Univariate Regression

import numpy as np
import scipy.odr.odrpack as odrpack
np.random.seed(1)

N = 100
x = np.linspace(0,10,N)
y = 3*x - 1 + np.random.random(N)
sx = np.random.random(N)
sy = np.random.random(N)

def f(B, x):
    return B*x + B
linear = odrpack.Model(f)
# mydata = odrpack.Data(x, y, wd=1./np.power(sx,2), we=1./np.power(sy,2))
mydata = odrpack.RealData(x, y, sx=sx, sy=sy)

myodr = odrpack.ODR(mydata, linear, beta0=)
myoutput = myodr.run()
myoutput.pprint()

Преобразование Fourier

Анализ Fourier помогает нам выразить функцию, как сумму периодических компонентов и восстановить сигнал из этих компонентов.

Давайте посмотрим на простой пример преобразования Fourier. Мы будем строить сумму двух синусов:

# Import Fast Fourier Transformation requirements
from scipy.fftpack import fft
import numpy as np

# Number of sample points
N = 600

# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)

# matplotlib for plotting purposes
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.grid()
plt.show()

SciPy Cluster

Clustering is the task of dividing the population or data points into a number of groups such that data points in the same groups are more similar to other data points in the same group and dissimilar to the data points in other groups. Each group which is formed from clustering is known as a cluster. There are two types of the cluster, which are:

  • Central
  • Hierarchy

Here we will see how to implement the K-means clustering algorithm which is one of the popular clustering algorithms. The k-means algorithm adjusts the classification of the observations into clusters and updates the cluster centroids until the position of the centroids is stable over successive iterations.

In the below implementation, we have used NumPy to generate two sets of random points. After joining both these sets, we whiten the data. Whitening normalizes the data and is an essential step before using k-means clustering. Finally, we use the kmeans functions and pass it the data and number of clustered we want.

import numpy as np
from scipy.cluster.vq import  kmeans, whiten
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("darkgrid")
# Create 50 datapoints in two clusters a and b
pts = 100
a = np.random.multivariate_normal(, 
                                  , ], 
                                  size=pts)
b = np.random.multivariate_normal(,
                                  , ],
                                  size=pts)
features = np.concatenate((a, b))
# Whiten data
whitened = whiten(features)
# Find 2 clusters in the data
codebook, distortion = kmeans(whitened, 2)
# Plot whitened data and cluster centers in red
plt.scatter(whitened, whitened)
plt.scatter(codebook, codebook, c='r')
plt.show()

Output:

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Adblock
detector