Similarly, an anti-Buddhabrot is the probability density distribution that a complex number does belong to a sequence generated from an element of the Mandelbrot set.

*An example of a Buddhabrot, on top, and of an anti-Buddhabrot, below, both obtained at once using 2000 iterations*

The Python code used is reported below.

```
# Computing, showing and saving raw data for Buddhabrot and anti-BBrot
# Dino Accoto 8/10 Oct 2021
# Credits: my starting point was the code published by prof. Susan Stepney,
# (retrieved on 8 Oct 2021) available here:
# http://susan-stepney.blogspot.com/2020/05/the-buddhabrot.html
# It needed some minor modifications
#to make it run on my Anaconda/Windows distribution.
# I added some image and file management functions to enable post-processing.
# With the following settings: Image_Size = 2048, Iterations = 2000,
# how_many_more = 20, the computation time on a PC with
# Intel Core i7-1065G7 CPU @ 1.30GHz - 1.50GHz, 16 Gb RAM
# is about 6 hours and 10 min.
# The code has been also tested on an Android smartphone with Pydroid 3.
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime # to create unique file names using timestamp
Image_Size = 2048 # Image Width and heigth (2048)
Iterations = 2000 # at least 500
resolution = 600 # image DPI (at least 300)
# The number of seed points is equal to the number
# of pixels (i.e. Image_Size**2) times 'how_many_more':
how_many_more = 20 # at least 10
def mandelbrot(c):
exp=2 # can play with exponent to get other shapes
z = c
sequence = [c]
for i in range(1, Iterations):
z = z ** exp + c
sequence += [z]
if abs(z) > 2: # the sequence escaped: can stop and record for Bud
return sequence, []
return [], sequence # sequence did not escape: recorded for ABud
def updateimage(img, seq):
for z in seq:
x, y = z.real, z.imag
ix, iy = int((x+2)/4*Image_Size), int((2-y)/4*Image_Size)
# Increment only if (ix,iy) is inside the image:
if (ix >= 0) and (ix < Image_Size) and (iy >= 0) and (iy < Image_Size):
img[ix,iy] += 1
def file_name(start_string):
filename = str(start_string)+'S'+str(Image_Size)+'I'+str(Iterations)+'R'+str(resolution)+'-'+str(datetime.now())
filename = filename.replace(':','-')
filename = filename.replace(' ','-')
filename = filename.replace('.','-')
return filename
def show_save(theimage,start_str): #show and save image
plt.axis('off')
plt.imshow(theimage, cmap='bone') #gist_earth cubehelix ocean plasma
plt.savefig(file_name(start_str), bbox_inches="tight", dpi=resolution)
plt.show()
################### Main Program ###########################
# Initialize the two images with 1s, so we can use log later, if needed
buddha = np.ones([Image_Size,Image_Size])
abuddha = np.ones([Image_Size,Image_Size])
for i in range(Image_Size*Image_Size*how_many_more):
# Option 1: seeds randomly picked in
# the 4x4 rectangle:
z = complex(np.random.uniform()*4-2, np.random.uniform()*4-2)
# Option 2: seeds randomly picked in a R=2 circle using polar coords
# (!but this option picks more seeds from near (0,0)):
# rho = np.random.uniform()*2
# theta = np.random.uniform()*2*np.pi
# z = complex(rho*np.cos(theta), rho*np.sin(theta))
# Option 3 to be used with Option 1: of all the seeds, just those in the
# R=2 circle are used. This highlights a circle around Bud/aBud:
# if (np.abs(z)<2):
traj, traja = mandelbrot(z)
updateimage(buddha,traj)
updateimage(abuddha,traja)
# Now the image has been created. Let's save it as a .npy binary file
# before further processing. This allows to experiment different cmaps
# and figure enhancing features.
# ==> To load the matrix from file: matr = np.load(file-name.pny)
np.save(file_name('B'),buddha)
np.save(file_name('A'),abuddha)
# Now some visual adjustments:
buddha = np.square(np.log(buddha))
abuddha = np.log(abuddha)
# Finally, let's show on screen and save on file as .png:
show_save(buddha,'B-')
show_save(abuddha,'a-')
```

]]>This ability, having an impact on the measured data, can alter the perception of the problem.

For example, a country with a low capability to detect infections could mistakenly appear, and consider itself, less affected by the epidemic than another country, perhaps less affected in absolute terms but with a better detection capacity.

We can expect that many heterogeneous factors affect the detection capability. Such factors may include: the organization and capillarity of health systems; the level of attention of the population to health-related issues; the spread of information; the ease of access to care.

It is difficult – perhaps impossible – to identify a single synthetic factor that can represent all these aspects in a synoptic way. Among the possible options,

We tried to evaluate whether there could be any link between the number of infections in European countries, expressed as the total number of cases per million inhabitants, and the life expectancy of that country.

A preliminary analysis showed that there is actually an exponential trend: the number of infections per million inhabitants grows exponentially with life expectancy in a given country.

The fitting was repeated over several days, from 18 to 28 March 2020, verifying a monotonous growth in the determination coefficient ().

The following figure shows the fitting with data related to 27 March 2020.

We found the following correlation law ():

(1)

where is life expectancy at birth and is the number of infections per million inhabitants.

To evaluate the predictive capacity of this equation, we proceeded as follows.

Since each Italian region is characterized by its own regional health system and specific demographic data, including life expectancy, the relationship (1) has been used to predict the number of infections.

The following figure compares predicted and real values (official data updated on 27/03/20).

As expected, the infection rate is lower than the real one only for the northern regions, most affected by the Covid-19 epidemic; the expected value is instead an underestimation of the real value for the central-southern regions, which probably benefited from the containment measures implemented by the Italian government, that slowed down the spread of the virus.

In conclusion, it is believed that there is a non-negligible correlation between the number of infections observed in a country and life expectancy in that country.

Furthermore, the relationship found has a certain predictive power, although the data are affected by a significant dispersion.

The origin of this relationship remains to be understood. It is probable that the effort for its understanding may shed light on the deeper implications of diversity in life expectancy among various countries.

]]>Epidemiologists and data scientists resort to multi-factorial models to predict the dynamics of the outbreak. Simpler tools can be used to look for trends just fitting epidemiological data.

The simplest models are based on either polynomial (second or higher order) or exponential functions. Both of them can provide a reasonable fitting during the first days of the outbreak, but not for long, as both of them diverge to infinity with time, while the total cases or deaths and infections should plateau to a value that cannot exceed the number of people in a given country.

A more accurate model is based on the ** Logistic function**, which belongs to the wider class of the so-called

With the Matlab script described hereafter, you will be able to generate fitting plots as the ones shown below, **using your own data**.

Let’s see how infection data can be fitted in **Matlab **using the Logistic function.

The analytical expression of our function is:

First of all, we need some data to fit.

**In this Excel file you find example data**, describing the Covid-19 outbreak in Italy in the past days. The structure of the file should be quite self-explaining.

If you want to use your own data, there are several reliable sources of information, such as the World Health Organization website, or this one, very convenient.

The Matlab script described in the following will show you how such data can be fitted using the Logistic function and its derivative.

```
% Matlab script to fit the number of total cases and infections using
% the logistic function: A/(1+exp(-B*(x-C))).
%
% Once the coefficients A, B e C are found, the (time) derivative is computed and the rate of
% change is compared with actual data about new cases and new deaths.
%
% Singapore, March 13 and 14, 2020
clear all
clc
close all
% "day" is the variable counting days in the observation interval.
% The range of "day" is from 1 to "Last_day"
Last_day = 60; % This number can be free
prev_day = 1:Last_day;
Data_from_spreadsheet = readtable("covid_data.xlsx");
% Let's check the size of the table read from Excel.
% "num_day" is the number of days considered.
[num_day,D2] = size(Data_from_spreadsheet);
% Let's convert the data in the table into a matrix "Data":
Data = Data_from_spreadsheet{1:num_day,2:5};
day = [1:num_day]';
today = max(day);
cases = Data(:,1);
new_cases = Data(:,2);
deaths = Data(:,3);
new_deaths = Data(:,4);
% =============================================================
% Let's find the coefficients of the logistic function by fitting data
% about total cases:
ft = fittype( 'a/(1+exp(-b*(x-c)))', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
% We need to provide some guess values to the algorithm to find
% the 3 coefficients (a, b, c).
% The fitting algorithm does not always converge. Choosing proper guess values
% is important.
opts.StartPoint = [70000 0.25 30];
% Data fitting:
[fitresult, gof] = fit( day, cases, ft, opts );
A = fitresult.a;
B = fitresult.b;
C = fitresult.c;
% Total cases according to fitting (logistic function):
prev = A./(1+exp(-B.*(prev_day-C)));
% New cases per day (derivative of logistic function, analytically computed):
derivative = A.*B.*exp(-B.*(prev_day-C))./(exp(-B.*(prev_day-C))+1).^2;
hold on
plot(prev_day,prev); %figure 1
plot(day,cases,'ro');
% Let's plot a horizontal line corresponding to the asymptot :
yline(A,'-.r');
% Let's report the R^2 on the plot:
text(today/3,max(prev)/2,strcat("R^2 = ",num2str(gof.rsquare)));
text(today/3,A*1.05,strcat("Total cases = ",num2str(round(A))));
%ylim([0 1.1*A]); %could be useful to set the range of values on the y-axis
title("Total cases (fitting with Logistic function)");
xlabel("Day");
hold off
figure
hold on
plot(prev_day,derivative,'b'); grid on %this creates figure 2
plot(day,new_cases,'ro');
hold off
xline(today,'-.r');
[Max,Index] = max(derivative); %this spots the max value
xline(Index);
text(Index,Max,strcat(" Peak: ",num2str(round(Max))," cases/day (day ",num2str(Index),")"));
%ylim([0 1.1*Max]);
title("New cases per day (derivative of Logistic function)");
xlabel("Day");
text(today,Max*0.1,strcat("Today: day ",num2str(today)));
% ======== let's repeat the same but for deaths ==========
ftM = fittype( 'a/(1+exp(-b*(x-c)))', 'independent', 'x', 'dependent', 'y' );
optsM = fitoptions( 'Method', 'NonlinearLeastSquares' );
optsM.Display = 'Off';
optsM.StartPoint = [10000 0.25 30];
% Data fitting:
[fitresultM, gofM] = fit( day, deaths, ftM, optsM );
AM = fitresultM.a;
BM = fitresultM.b;
CM = fitresultM.c;
% Total deaths according to fitting
prevM = AM./(1+exp(-BM.*(prev_day-CM)));
derivativeM = AM.*BM.*exp(-BM.*(prev_day-CM))./(exp(-BM.*(prev_day-CM))+1).^2;
figure
hold on
plot(prev_day,prevM); % this creates fig 3
plot(day,deaths,'ro');
% Horizontal asymptote:
yline(AM,'-.r');
text(today/3,max(prevM)/2,strcat("R^2 = ",num2str(gofM.rsquare)));
text(today/3,AM*1.05,strcat("Total deaths = ",num2str(round(AM))));
%ylim([0 1.1*AM]);
title("Total deaths (fitting with Logistic function)");
xlabel("Day")
hold off
figure
hold on
plot(prev_day,derivativeM,'b'); grid on % this creates fig 4
plot(day,new_deaths,'ro');
hold off
xline(today,'-.r');
[MaxM,IndexM] = max(derivativeM);
xline(IndexM);
text(IndexM,MaxM,strcat(" Peak: ",num2str(round(MaxM))," deaths/day (day ",num2str(IndexM),")"));
%ylim([0 1.1*MaxM]);
title("New deaths/day (derivative of Logistic function)");
xlabel("Day");
text(today,MaxM*0.1,strcat("Today: day ",num2str(today)));
```

]]>