Стабилизация продуктовой метрики Удержание (Retention)
Новый оператор тензорного произведения и анализ нелинейных систем с полиномиальными функциями пространства состояний
(см. Минитаева А.М. Новый оператор тензорного произведения и анализ нелинейных систем с полиномиальными функциями пространства состояний // Информационно-технологический вестник. 2024. № 3.)
Новый оператор умножения многоиндексных объектов
Оператор производит умножение с векторным результатом по одному измерению тензоров и фробениусово скалярное умножение по остальным измерениям. Такой бинарный оператор тензорной алгебры можно определить следующим образом:
\begin{equation}\label{eq19}
\begin{gathered}
{A^{d + 1}}\mathop \otimes \limits^\mu {B^{d + 1}} = \
= ({A_1},{A_2},...,{A_d})\mathop \otimes \limits^\mu ({B_1},{B_2},...,{B_d}) = \
= \left( {\begin{array}{*{20}{c}}
{{A_1}\mathop \times \limits^{\varphi \rho } {B_1}} \
{{A_2}\mathop \times \limits^{\varphi \rho } {B_2}} \
{...} \
{{A_d}\mathop \times \limits^{\varphi \rho } {B_d}}
\end{array}} \right), \hfill \
{ A,B} \in { {\mathbb{R}{{d{d + 1}}}},{{\mathbb{X}}{{d{d + 1}}}}} ,
\end{gathered}
\end
(1)
где - это тензоры ранга ; , - многоиндексные объекты ранга .
Оператор по выражению (1) можно реализовать в алгоритмических языках программирования, например, Python с подключением библиотеки numpy.
В качестве координаты на пространстве можно использовать булевский тип данных, т.к. упорядоченный набор индексов тензора на пространстве задает (индексы в определенном измерении являются степенями соответствующей этому измерению переменной) одночлен переменных с максимальной степенью из вектора размерности .
Ниже представлен код Python, выполняющий умножение многоиндексных объектов на основе данного определения:
import numpy as np
def frobenius_inner_product(A, B):
# Calculate the Frobenius inner product.
return np.trace(np.dot(np.transpose(A), B))
def mu_tensor_product(A, B):
# Calculate the novel tensor product.
# Get the dimensions of A and B
shape_A = A.shape
shape_B = B.shape
# Ensure A and B have the same shape
if shape_A != shape_B:
raise ValueError("Shapes of A and B must match")
# Initialize the result vector
result = np.zeros(shape_A[0])
# Perform the tensor product operation
for i in range(shape_A[0]):
result[i] = frobenius_inner_product(A[i], B[i])
return result
# Example usage
if __name__ == "__main__":
# Define two example multi-index objects A and B
A = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]])
B = np.array([[[.9,.8], [.7,.6]], [[.5,.4], [.3,.2]]])
# Perform the mu_tensor_product
result = mu_tensor_product(A, B)
# Print the result
print("Result:")
print(result)
# [7. 8.6]
Тензорная форма записи -го уравнения состояния нелинейной системы управления с полиномиальными функциями пространства состояний
Полином содержит членов:
[t,x = t,{x_1},{x_2} = {a_{i,2,0}}x_1^2 + {a_{i,1,1}}{x_1}{x_2} + {a_{i,0,2}}x_2^2 + {a_{i,1,0}}{x_1} + {a_{i,0,1}}{x_2}.]
Назовем
\begin{equation}\label{eq8}
{X_C} = \left( {\begin{array}{*{20}{c}}
0&{{x_2}}&{x_2^2} \
{{x_1}}&{{x_1}{x_2}}&0 \
{x_1^2}&0&0
\end{array}} \right)
\end
матрицей сочетаний переменных состояния, а
\begin{equation}\label{eq9}
{A_i} = \left( {\begin{array}{*{20}{c}}
0&{{a_{i,0,1}}}&{{a_{i,0,2}}} \
{{a_{i,1,0}}}&{{a_{i,1,1}}}&0 \
{{a_{i,2,0}}}&0&0
\end{array}} \right)
\end
матрицей коэффициентов полинома состояния -го уравнения состояния.
Для всей системы из уравнений можно говорить о многоиндексном объекте ранга 3:
\begin{equation}\label{eq10}
A = \left( {\left( {\begin{array}{{20}{c}}
0&{{a_{1,0,1}}}&{{a_{1,0,2}}} \
{{a_{1,1,0}}}&{{a_{1,1,1}}}&0 \
{{a_{1,2,0}}}&0&0
\end{array}} \right)\left( {\begin{array}{{c}}
0&{{a_{2,0,1}}}&{{a_{2,0,2}}} \
{{a_{2,1,0}}}&{{a_{2,1,1}}}&0 \
{{a_{2,2,0}}}&0&0
\end{array}} \right)} \right),
\end{equation}
где .
Также создадим многоиндексный объект ранга 3 из матрицы , продублировав ее в качестве его сечения:
\begin{equation}\label{eq11}
X = \left( {\left( {\begin{array}{{20}{c}}
0&{{x_2}}&{x_2^2} \
{{x_1}}&{{x_1}{x_2}}&0 \
{x_1^2}&0&0
\end{array}} \right),\left( {\begin{array}{{c}}
0&{{x_2}}&{x_2^2} \
{{x_1}}&{{x_1}{x_2}}&0 \
{x_1^2}&0&0
\end{array}} \right)} \right),
\end
где ,
- это бесконечное множество одночленов для переменных вектора .
Теперь можно записать уравнение САУ в форме многоиндексных объектов:
\begin{equation}\label{eq12}
{\Delta x_i}[n + 1] = {A_i}\mathop \times \limits^{\varphi \rho } {X_C} + {B_i}u[n],
\end
где - бинарный оператор фробениусова скалярного произведения (или стандартного скалярного произведения), определенного на пространствах любого ранга , как сумма произведений элементов, чьи индексы полностью совпадают: . Фробениусово скалярное произведение может также быть вычислено как след (сумма элементов главной диагонали) произведения транспонированного тензора и тензора .
Форма записи с использованием многоиндексных объектов системы уравнений пространства состояний нелинейной системы управления с дискретным временем
Обозначим:
-
- многоиндексный объект сочетания переменных состояния и управления ранга , где - количество переменных состояния, и - количество переменных управления.
-
- многоиндексный объект (ранга ) коэффициентов полинома состояния -го уравнения состояния системы.
-
- многоиндексный объект, получаемый дублированием в качестве своего сечения, ранга .
-
- многоиндексный объект (ранга ) коэффициентов полиномиальной системы уравнений нелинейной системы управления.
Теперь все по уравнения пространства состояний нелинейной системы управления с дискретным временем можем записать в форме многоиндексных объектов::
\begin{equation}\label{eq22}\Delta
x[n + 1] = A\mathop \otimes \limits^\mu W.
\end
Уравнение нелинейной системы управления еще компактнее векторной формы уравнения линейной системы и зачастую может быть удобнее для представления данных в ЭВМ и алгоритмических языках программирования, а также реализует новый тензорный подход для моделирования нелинейных систем управления с дискретным временем.
Моделирование САУ
Управляющий сигнал и модель пространства состояний
# Управляющий сигнал (Ступенька)
u_max = 0.0
u_time = 1.0
A = np.array([[[0, 0, 1],
[0, 0, 0],
[0, 0, 0]],
[[0, 0, 0],
[0, 0, 0],
[1, 0, 0]]])
t=0
x1=1
x2=2
X = np.array([[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]],
[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]]])
# Perform the mu_tensor_product
result = -mu_tensor_product(A, X)
# Print the result
print("Result:")
print(result)
# Создаем списки значений для обеих функций
xs1 = []
ys1 = []
ts1 = []
fs1 = []
cs1 = []
xs2 = []
ys2 = []
ts2 = []
fs2 = []
cs2 = []
ts3 = []
xs1.append(t)
ys1.append(result[1])
ts1.append(x1)
ts3.append(x1+1)
fs1.append((x1+1)/2)
xs2.append(t)
ys2.append(result[0])
ts2.append(x2)
fs2.append(x2/2)
cs1.append(0)
cs2.append(0)
step=0.01
sum_1 = 0
sum_2 = 0
sum_d1 = 0
sum_d2 = 0
sum_c1 = 0
sum_c2 = 0
sum_1 += ((x1+1)/2)*step
sum_2 += (x2/2)*step
sum_d1 += result[1]*step
sum_d2 += result[0]*step
sum_c1 += 0
sum_c2 += 0
x1 += result[1]*step
x2 += result[0]*step
X = np.array([[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]],
[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]]])
# Perform the mu_tensor_product
result = -mu_tensor_product(A, X)
xs1.append(t+step)
ys1.append(result[1])
ts1.append(x1)
ts3.append(x1+1)
fs1.append((x1+1)/2)
xs2.append(t+step)
ys2.append(result[0])
ts2.append(x2)
fs2.append(x2/2)
cs1.append(0)
cs2.append(0)
sum_1 += ((x1+1)/2)*step
sum_2 += (x2/2)*step
sum_d1 += result[1]*step
sum_d2 += result[0]*step
sum_c1 += 0
sum_c2 += 0
# Print the result
print("Result:")
print(result)
Симуляция
import matplotlib.pyplot as plt
# Границы изменения переменной t
t_min = 0.0
t_max = 6.0
# Шаг для вычислений
step = 0.01
# Период вычисления оттока
churn_period = 0.23 # 1 неделя: 0.23=7*12/365.25
# Управляющий сигнал
u = 0.0
t = t_min + step
while t<=t_max:
if t >= u_time:
u = u_max
x1 += result[1]*step
x2 += result[0]*step
X = np.array([[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]],
[[0, x2, x2 ** 2],
[x1, x1 * x2, 0],
[x1 ** 2, 0, 0]]])
# Perform the mu_tensor_product
result = -mu_tensor_product(A, X) + [u, u]
xs1.append(t)
ys1.append(result[1])
ts1.append(x1)
ts3.append(x1+1)
fs1.append((x1+1)/2)
xs2.append(t)
ys2.append(result[0])
ts2.append(x2)
fs2.append(x2/2)
if t >= churn_period:
prev = int((t-churn_period)/step)
cs1.append((fs1[prev]-(x1+1)/2)/fs1[prev])
cs2.append((fs2[prev]-x2/2)/fs2[prev])
sum_c1 += step*(fs1[prev]-(x1+1)/2)/fs1[prev]
sum_c2 += step*(fs2[prev]-x2/2)/fs2[prev]
else:
cs1.append(0)
cs2.append(0)
sum_1 += ((x1+1)/2)*step
sum_2 += (x2/2)*step
sum_d1 += result[1]*step
sum_d2 += result[0]*step
t+=step
# Построение графика
plt.figure()
plt.title('Графики при н.у. x=1, y=2')
plt.grid(True)
plt.plot(xs1, ts1, label='x')
plt.plot(xs2, ts2, label='y')
plt.plot(xs1, ts3, linestyle='dashed', label='x+1')
plt.legend()
plt.xlabel('t')
plt.ylabel('x, y')
plt.show()
plt.figure()
plt.title('Графики кривых удержания')
plt.grid(True)
plt.plot(xs1, fs1, label='x1')
plt.plot(xs2, fs2, label='x2')
plt.legend()
plt.xlabel('t')
plt.ylabel('Retention')
plt.show()
print('Площадь под графиком x1: ',sum_1)
print('Площадь под графиком x2: ',sum_2)
plt.figure()
plt.title('Графики производных')
plt.grid(True)
plt.plot(xs1, ys1, label='x\'')
plt.plot(xs2, ys2, label='y\'')
plt.legend()
plt.xlabel('t')
plt.ylabel('x\', y\'')
plt.show()
print("Значения производных в точке t= "+str(u_time))
print("x\'("+str(u_time)+") = "+"%.2f" % ys1[int(u_time/step)-1])
print("y\'("+str(u_time)+") = "+"%.2f" % ys2[int(u_time/step)-1])
print()
print('Площадь под графиком x\': ',sum_d1)
print('Площадь под графиком y\': ',sum_d2)
ut=[0,u_time,u_time,t_max]
us=[0,0,u_max,u_max]
plt.figure()
plt.title('Управляющий сигнал')
plt.grid(True)
plt.plot(ut, us, label='u')
plt.legend()
plt.xlabel('t')
plt.ylabel('Ступенька')
plt.show()
del cs1[0:int(churn_period/step)+1]
del cs2[0:int(churn_period/step)+1]
for i in range(int(churn_period/step)+1):
cs1.append(cs1[int((t_max-churn_period)/step)])
cs2.append(cs2[int((t_max-churn_period)/step)])
plt.figure()
plt.title('Графики оттока в '+str(churn_period)+' ед.врем.')
plt.grid(True)
plt.plot(xs1, cs1, label='Churn 1')
plt.plot(xs2, cs2, label='Churn 2')
plt.legend()
plt.xlabel('t')
plt.ylabel('Churn')
plt.show()
print('Площадь под графиком Churn 1: ',sum_c1)
print('Площадь под графиком Churn 2: ',sum_c2)
.png)
.png)
.png)
.png)
Графики
plt.figure()
plt.title('Фазовый портрет процессов удержания x1, x2')
plt.grid(True)
plt.plot(fs1, fs2, label='x2(x1)')
plt.legend()
plt.xlabel('x1')
plt.ylabel('x2')
plt.arrow (x= 0.9 , y= 0.8 , dx= -0.05 , dy= -0.1 , width= .01 )
plt.show()
plt.figure()
plt.title('Фазовые портреты производной и процесса удержания x\', x')
plt.grid(True)
plt.plot(fs1, ys1, label='x1\'(x1)')
plt.plot(fs2, ys2, label='x2\'(x2)')
plt.legend()
plt.xlabel('x')
plt.ylabel('x\'')
plt.arrow (x= 0.9 , y= -2.0 , dx= -0.1 , dy= 0.4 , width= .02 )
plt.show()
plt.figure()
plt.title('Фазовые портреты оттока и удержания')
plt.grid(True)
plt.plot(fs1, cs1, label='Churn1(x1)')
plt.plot(fs2, cs2, label='Churn2(x2)')
plt.legend()
plt.xlabel('Retention')
plt.ylabel('Churn')
plt.arrow (x= 0.9 , y= 0.2 , dx= -0.1 , dy= -0.03 , width= .01 )
plt.show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.title('Фазовые портреты оттоков и удержаний')
plt.grid(True)
cs=np.ones((len(cs1), 10))
c=np.ones(10)
fs_1=np.ones((len(fs1), 10))
fs_2=np.ones((len(fs2), 10))
fc_1=np.ones(10)
fc_2=np.ones(10)
for i in range(len(cs1)):
c = np.mgrid[cs1[i]:cs2[i]:10j]
fc_1 = np.mgrid[fs1[i]:fs1[i]:10j]
fc_2 = np.mgrid[fs2[i]:fs2[i]:10j]
cs[i] = c
fs_1[i] = fc_1
fs_2[i] = fc_2
ax.plot_surface(fs_1, fs_2, cs, cmap='plasma')
ax.plot(fs1, fs2, cs1, label='Churn1(x1, x2)')
ax.plot(fs1, fs2, cs2, label='Churn2(x2, x2)')
plt.legend()
plt.xlabel('Retention-1')
plt.ylabel('Retention-2')
plt.show()
.png)
.png)
.png)
Фазовые портреты
import random
import math
import pandas
import pandas as pd
def pdf(k: int, rnd_list: list) -> pandas.DataFrame:
"""
Получает кривую плотности распределения вероятности
:param k: количество интервалов разбиения гистограммы
:param rnd_list: случайный процесс
:return: pandas.DataFrame
"""
pdf_x = [] # Координаты по оси абсцисс
pdf_y = [] # Координаты по оси ординат
n = len(rnd_list) # количество элементов в рассматриваемой выборке
h = (max(rnd_list) - min(rnd_list)) / k # ширина одного интервала
a = min(rnd_list) # минимальное значение в рассматриваемой выборке
for i in range(0, k): # проход по интервалам
count = 0
for j in rnd_list: # подсчет количества вхождений значений из выборки в данный интервал
if (a + i * h) < j < (a + (i * h) + h):
count = count + 1
pdf_x.append(a + i * h + h / 2) # координата по оси абсцисс полученной кривой плотности распределения
# вероятности
pdf_y.append(count / (n * h)) # координата по оси ординат полученной кривой плотности распределения
# вероятности
d = {'x': pdf_x, 'y': pdf_y}
return pd.DataFrame(d)
rrand = fs1 # числа с распределением удержания x1
r_pdf = pdf(100, rrand) # строим плотность распределения вероятности по полученной выборке
grand = fs2 # числа с распределением удержания x2
g_pdf = pdf(100, grand) # строим плотность распределения вероятности по полученной выборке
wrand = cs1 # числа с распределением оттока Churn1
w_pdf = pdf(100, wrand) # строим плотность распределения вероятности по полученной выборке
exprand = cs2 # числа с распределением оттока Churn2
e_pdf = pdf(100, exprand) # строим плотность распределения вероятности по полученной выборке
# Построение графиков
fig, ax1 = plt.subplots(nrows=2, ncols=2) # Создаем фигуру и четыре системы координат на фигуре
ax1[0, 0].plot(r_pdf['x'], r_pdf['y'], 'b-', label='Retention-1') # Своя оценка плотности распределения
ax1[0, 0].set_ylabel('PDF(x)') # Ось ординат
ax1[0, 0].set_xlim(xmin=0.5, xmax=1) # Минимальное и максимальное значение по оси Х
ax1[0, 0].legend() # Отображаем легенду для данной системы координат
ax1[0, 0].grid() # Отображаем сетку на системе координат
ax1[0, 1].plot(g_pdf['x'], g_pdf['y'], 'r-', label='Retention-2') # Своя оценка плотности распределения
ax1[0, 1].set_xlim(xmin=0, xmax=1) # Минимальное и максимальное значение по оси Х
ax1[0, 1].legend() # Отображаем легенду для данной системы координат
ax1[0, 1].grid() # включение отображение сетки
ax1[1, 0].plot(w_pdf['x'], w_pdf['y'], 'b-', label='Churn-1')
ax1[1, 0].set_ylabel('PDF(churn)') # ось ординат
ax1[1, 0].set_xlim(xmin=-0.05, xmax=0.1)
ax1[1, 0].legend()
ax1[1, 0].grid() # включение отображение сетки
ax1[1, 1].plot(e_pdf['x'], e_pdf['y'], 'r-', label='Churn-2')
ax1[1, 1].set_xlim(xmin=0, xmax=0.35)
ax1[1, 1].legend()
ax1[1, 1].grid() # включение отображение сетки
Плотность распределения вероятности
Имитационное моделирование удержания с помощью марковского процесса с непрерывным временем
Поток событий — это последовательность однородных событий, наступающих одно за другим в случайные промежутки времени.
Интенсивность потока λ — это среднее число событий в единицу времени. Интенсивность потока можно рассчитать экспериментально по формуле: λ = N/Tн, где N — число событий, произошедших за время наблюдения Tн.
За эталон потока в моделировании принято брать пуассоновский поток.
Пуассоновский поток — это ординарный поток без последействия.
Вероятность того, что за интервал времени (t0, t0 + τ) произойдет m событий, определяется из закона Пуассона
где a — параметр Пуассона.
Если λ(t) = const(t), то это стационарный поток Пуассона (простейший). В этом случае a = λ · t.
Если λ = var(t), то это нестационарный поток Пуассона.
Рисунок иллюстрирует зависимость P0 от времени. Чем больше значение λ, тем круче идет график, то есть быстрее убывает вероятность.
В нашем случае для двух процессов удержания и показатель λ равен минус значению производной в начальной точке, т.е. и .
Тогда случайный интервал между двумя событиями в этих двух потоках и , где R - это равномерно распределенное от 0 до 1 случайное число, каждый раз новое.
Пусть в начальный момент времени мы имеем 1000 клиентов (500 пришедших в нулевой момент времени и 500 в ядре аудитории) в процессе удержания и 1000 клиентов в процессе удержания . По формуле выше смоделируем моменты времени оттока каждого клиента.
import random
import math
import matplotlib.pyplot as plt
l1 = 1
l2 = 4
N1 = 500
N2 = 1000
t1 = []
t2 = []
for i in range(N1):
t1.append(-math.log(random.random())/l1)
for i in range(N2):
t2.append(-math.log(random.random())/l2)
dt = 0.01
a = 0
b = 2
t = a
ts = []
counts1 = []
counts2 = []
for i in range(int(b/dt+1)):
cnt1 = 0
cnt2 = 0
ts.append(t)
for j in t1:
if j > t:
cnt1+=1
for k in t2:
if k > t:
cnt2+=1
counts1.append(cnt1+N1)
counts2.append(cnt2)
t += dt
plt.figure()
plt.title('Графики кривых удержания (имитационно)')
plt.grid(True)
plt.plot(ts, counts1, label='x1')
plt.plot(ts, counts2, label='x2')
plt.legend()
plt.xlabel('t')
plt.ylabel('Retention')
plt.show()
plt.figure()
plt.title('Графики кривых удержания (аналитически)')
plt.grid(True)
plt.plot(xs1, fs1, label='x1')
plt.plot(xs2, fs2, label='x2')
plt.legend()
plt.xlabel('t')
plt.ylabel('Retention')
plt.show()
.png)
Разница во временных границах имитационных и аналитических графиков объясняется тем, что мы использовали натуральный логарифм в формулах и для имитационного моделирования марковских процессов.