---
title: "Λογιστική Αύξηση Πληθυσμού"
subtitle: "Όταν ο χώρος αρχίζει να μετράει: από την εκθετική στη λογιστική εξίσωση"
author: "Εργαστήριο Πληθυσμιακής Οικολογίας — Τμήμα Βιολογίας"
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    number_sections: true
    theme: flatly
    highlight: tango
    code_folding: show
    df_print: paged
  pdf_document:
    toc: true
    number_sections: true
    latex_engine: xelatex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 7,
  fig.height = 4.5,
  fig.align = "center"
)
library(ggplot2)
library(deSolve)
theme_set(theme_minimal(base_size = 12))
```

# Πού μείναμε

Στο προηγούμενο εργαστήριο καταλήξαμε στον τύπο της εκθετικής αύξησης:

$$
\frac{dN}{dt} = r N \quad\Longrightarrow\quad N(t) = N_0\, e^{rt}
$$

Με $N_0 = 100$ νούφαρα και $r = 0{,}1$/ημέρα, ας δούμε τι προβλέπει το μοντέλο σε βάθος χρόνου:

```{r exp-projection}
N0 <- 100
r  <- 0.1

cat("Σε 100 ημέρες:", format(N0 * exp(r * 100),  big.mark = "."), "νούφαρα\n")
cat("Σε 200 ημέρες:", format(N0 * exp(r * 200),  big.mark = "."), "νούφαρα\n")
cat("Σε 365 ημέρες:", format(N0 * exp(r * 365),  big.mark = "."), "νούφαρα\n")
```

**Είναι παράλογο.** Σε λιγότερο από έναν χρόνο, η μικρή μας λίμνη υποτίθεται ότι θα έχει περισσότερα νούφαρα από όσα κύτταρα περιέχει το ανθρώπινο σώμα.

Το πρόβλημα είναι ότι το εκθετικό μοντέλο **αγνοεί τον χώρο, το φως, τα θρεπτικά συστατικά και τον ανταγωνισμό**. Σε αυτό το εργαστήριο θα διορθώσουμε αυτή την έλλειψη.

---

# Η ιδέα της φέρουσας ικανότητας $K$

Κάθε πραγματικό περιβάλλον έχει ένα **όριο** στον αριθμό ατόμων που μπορεί να συντηρήσει. Ονομάζουμε αυτό το όριο **φέρουσα ικανότητα** (carrying capacity) και το συμβολίζουμε με $K$.

Για τη λίμνη μας, ας υποθέσουμε:

```{r set-K}
K <- 10000   # μέγιστος αριθμός νούφαρων που χωρά η λίμνη
```

> *Σημείωση επιλογής:* Η τιμή $K = 10{.}000$ είναι ενδεικτική και επιλέχθηκε για λόγους σαφήνειας στα διαγράμματα. Σε πραγματική μελέτη το $K$ εκτιμάται από δεδομένα ή από τα φυσικά χαρακτηριστικά του οικοτόπου (επιφάνεια λίμνης, μέσο μέγεθος νούφαρου, διαθέσιμο φως κλπ).

## Το $K$ δεν είναι απλώς "ταβάνι"

Η αφελής ιδέα θα ήταν: *"αφήνουμε τον πληθυσμό να αυξάνεται εκθετικά μέχρι το $K$, και μετά τον σταματάμε."* Αυτό όμως δεν συμβαίνει στη φύση — η αύξηση **επιβραδύνεται σταδιακά** καθώς πλησιάζουμε στο $K$, *πολύ πριν* το φτάσουμε.

Γιατί;

---

# Ένας μηχανισμός: γιατί το $K$ μειώνει τον ρυθμό

Φανταστείτε τη λίμνη μας με τα νούφαρα. Όταν ο πληθυσμός είναι μικρός, κάθε νούφαρο έχει άπλετο χώρο γύρω του και μπορεί να αναπαραχθεί ελεύθερα. Όμως καθώς ο πληθυσμός μεγαλώνει, **τα νούφαρα στο εσωτερικό της συστάδας περιβάλλονται από άλλα νούφαρα** και δυσκολεύονται να επεκταθούν. Μόνο τα νούφαρα στην *περιφέρεια* έχουν χώρο. Όσο πιο γεμάτη η λίμνη, τόσο μικρότερο το ποσοστό των νούφαρων που μπορούν πραγματικά να αναπαραχθούν.

Αυτό είναι **ένα παράδειγμα** μηχανισμού. Στην πραγματικότητα μπορούν να συντρέχουν πολλοί λόγοι:

- **Χώρος / φως:** τα νούφαρα στο εσωτερικό σκιάζονται.
- **Θρεπτικά:** ο ανταγωνισμός για ριζικά θρεπτικά συστατικά εντείνεται.
- **Παθογόνα:** σε πυκνές συστάδες οι ασθένειες εξαπλώνονται γρήγορα.
- **Τοξικά παραπροϊόντα:** συσσωρεύονται όταν ο πληθυσμός πυκνώνει.

Όλα αυτά μοιράζονται το ίδιο **ποιοτικό** χαρακτηριστικό: ο **κατά κεφαλήν ρυθμός αναπαραγωγής μειώνεται** καθώς ο πληθυσμός πλησιάζει το $K$.

> Άρα το $K$ δεν είναι απλώς ένα όριο που χτυπάμε επάνω. Το $K$ **τροποποιεί τον ίδιο τον ρυθμό μεταβολής**.

---

# Από τη φυσική διαίσθηση στη μαθηματική μορφή

## Τι θέλουμε από τον "πραγματοποιημένο" κατά κεφαλήν ρυθμό

Ορίζουμε τον **πραγματοποιημένο** (realized) κατά κεφαλήν ρυθμό αύξησης ως μια συνάρτηση $r_{\text{eff}}(N)$ του τρέχοντος πληθυσμού. Δύο φυσικές απαιτήσεις:

1. Όταν $N$ είναι **μικρό**, υπάρχει άπλετος χώρος → $r_{\text{eff}}(N) \to r$ *(ο εγγενής ρυθμός)*.
2. Όταν $N \to K$, ο χώρος εξαντλείται → $r_{\text{eff}}(N) \to 0$.

## Η απλούστερη συνάρτηση που ικανοποιεί αυτά: γραμμική

Η πιο απλή συνάρτηση που πιάνει τα δύο άκρα είναι μια **ευθεία γραμμή** που πέφτει από $r$ (στο $N=0$) στο $0$ (στο $N=K$):

$$
\boxed{\;r_{\text{eff}}(N) \;=\; r\left(1 - \frac{N}{K}\right)\;}
$$

> **Σημείωση επιλογής μοντέλου (σημαντική):** Η γραμμική μορφή είναι μια *μοντελοποιητική παραδοχή*, όχι λογική αναγκαιότητα. Είναι η απλούστερη συνάρτηση που πιάνει τη σωστή ποιοτική συμπεριφορά. Άλλες δυνατές μορφές (θα τις δούμε στο τέλος):
>
> - **θ-λογιστική:** $r_{\text{eff}}(N) = r\bigl(1 - (N/K)^{\theta}\bigr)$
> - **Gompertz:** $r_{\text{eff}}(N) = r\,\ln(K/N)$
>
> Όλες δίνουν παρόμοια σιγμοειδή συμπεριφορά, αλλά με διαφορετική γεωμετρία.

```{r per-capita-comparison}
N_grid <- seq(0, K, length.out = 200)
df_pc <- data.frame(
  N        = rep(N_grid, 2),
  rate     = c(rep(r, length(N_grid)),                # σταθερό r (εκθετικό)
               r * (1 - N_grid/K)),                   # γραμμικά μειούμενο (λογιστικό)
  Μοντέλο  = factor(rep(c("Εκθετικό:  r σταθερό",
                          "Λογιστικό: r·(1-N/K)"),
                        each = length(N_grid)))
)

ggplot(df_pc, aes(x = N, y = rate, colour = Μοντέλο)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(
    title = "Κατά κεφαλήν ρυθμός αύξησης ως συνάρτηση του πληθυσμού",
    subtitle = "Η ουσιαστική διαφορά μεταξύ εκθετικού και λογιστικού μοντέλου",
    x = "Πληθυσμός N",
    y = expression(r[eff](N))
  ) +
  theme(legend.position = "bottom")
```

---

# Η λογιστική εξίσωση

Ο συνολικός ρυθμός μεταβολής του πληθυσμού είναι:

$$
\frac{dN}{dt} \;=\; r_{\text{eff}}(N) \cdot N \;=\; r\,N\left(1 - \frac{N}{K}\right)
$$

Αυτή είναι η **λογιστική εξίσωση** (Verhulst, 1838).

$$
\boxed{\;\frac{dN}{dt} \;=\; r\,N\left(1 - \frac{N}{K}\right)\;}
$$

---

# Ποιοτική ανάλυση *πριν* λύσουμε

Ένα από τα πιο χρήσιμα εργαλεία στα δυναμικά συστήματα είναι να εξετάζουμε τη συμπεριφορά **πριν** λύσουμε αναλυτικά. Πολύ συχνά μάλιστα η λύση δεν υπάρχει αναλυτικά — αλλά η ποιοτική εικόνα μάς αρκεί.

## Σημεία ισορροπίας

Όπου $\dfrac{dN}{dt} = 0$ ο πληθυσμός μένει στάσιμος. Λύνουμε:

$$
r\,N\left(1 - \frac{N}{K}\right) = 0 \quad\Longrightarrow\quad N^* = 0 \;\text{ ή }\; N^* = K
$$

- $N^* = 0$: άδεια λίμνη (πληθυσμός εξαφανισμένος).
- $N^* = K$: γεμάτη λίμνη (κορεσμός).

## Ευστάθεια σημείων ισορροπίας

- Αν ξεκινήσουμε **κοντά** στο $0$ (π.χ. $N=1$): $\frac{dN}{dt} > 0$, ο πληθυσμός *μεγαλώνει* — απομακρύνεται από το 0. Άρα $N^*=0$ είναι **ασταθές**.
- Αν ξεκινήσουμε **κοντά** στο $K$ (π.χ. $N=0{,}9K$): $\frac{dN}{dt} > 0$, ο πληθυσμός *πλησιάζει* το $K$.
- Αν ξεκινήσουμε **πάνω** από $K$ (π.χ. $N=1{,}5K$): $\frac{dN}{dt} < 0$, ο πληθυσμός *μειώνεται* προς το $K$.

Άρα $N^*=K$ είναι **ευσταθές** — είναι ένας *ελκυστής* του συστήματος.

## Φασικό διάγραμμα: $\dfrac{dN}{dt}$ ως προς $N$

```{r phase-portrait}
N_grid <- seq(-500, K * 1.4, length.out = 300)
dN_dt  <- r * N_grid * (1 - N_grid/K)
df_phase <- data.frame(N = N_grid, dN = dN_dt)

# Σημεία ισορροπίας
eq_unstable <- data.frame(N = 0, dN = 0)
eq_stable   <- data.frame(N = K, dN = 0)

# Βέλη ροής στον άξονα N: το πρόσημο του (xend - x) δίνει την κατεύθυνση κίνησης
arrows_df <- data.frame(
  x    = c(K/4,        3*K/4,        1.2*K),
  xend = c(K/4 + 800,  3*K/4 + 800,  1.2*K - 800),
  y    = 0, yend = 0
)

ggplot(df_phase, aes(x = N, y = dN)) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(linewidth = 1, colour = "darkgreen") +
  geom_point(data = eq_unstable, aes(x = N, y = dN),
             shape = 21, fill = "white", colour = "darkgreen",
             size = 4, stroke = 1.2) +
  geom_point(data = eq_stable, aes(x = N, y = dN),
             colour = "darkgreen", size = 4) +
  geom_segment(data = arrows_df,
               aes(x = x, xend = xend, y = y, yend = yend),
               arrow = arrow(length = unit(0.25, "cm")),
               colour = "tomato", linewidth = 0.9,
               inherit.aes = FALSE) +
  annotate("text", x = 0,    y = -150, label = "N* = 0\n(ασταθές)",  vjust = 1) +
  annotate("text", x = K,    y = -150, label = "N* = K\n(ευσταθές)", vjust = 1) +
  annotate("text", x = K/2,  y = max(dN_dt) + 30,
           label = "μέγιστος ρυθμός\nστο N = K/2", vjust = 0) +
  labs(
    title = "Φασικό διάγραμμα της λογιστικής εξίσωσης",
    subtitle = "Τα βέλη δείχνουν την κατεύθυνση μεταβολής του N",
    x = "N", y = "dN/dt = r"
  )
```

Η παραβολή αναπαριστά **όλη τη δυναμική του συστήματος** σε μια εικόνα:

- Όπου η καμπύλη είναι **πάνω** από τον άξονα → $\frac{dN}{dt} > 0$ → ο πληθυσμός αυξάνεται.
- Όπου η καμπύλη είναι **κάτω** από τον άξονα → ο πληθυσμός μειώνεται.
- Στις *ρίζες* έχουμε σημεία ισορροπίας.
- Η **κορυφή** της παραβολής (στο $N = K/2$) είναι ο μέγιστος ρυθμός αύξησης — εκεί ο πληθυσμός μεγαλώνει πιο γρήγορα. Αυτό είναι και το **σημείο καμπής** της χρονικής τροχιάς.

---

# Αναλυτική λύση *(προαιρετικό για μαθηματικά διατεθειμένους)*

Για όσους θέλουν να δουν τη λύση: ξεκινάμε με χωρισμό μεταβλητών:

$$
\frac{dN}{N\,(1 - N/K)} = r\,dt
\quad\Longleftrightarrow\quad
\frac{K\,dN}{N\,(K - N)} = r\,dt
$$

**Παραγοντική ανάπτυξη** (partial fractions):

$$
\frac{K}{N(K-N)} = \frac{1}{N} + \frac{1}{K-N}
$$

Άρα:

$$
\int\!\left(\frac{1}{N} + \frac{1}{K-N}\right) dN = \int r\,dt
\;\Longleftrightarrow\;
\ln N - \ln(K-N) = rt + C
$$

$$
\ln\!\left(\frac{N}{K-N}\right) = rt + C
\;\Longrightarrow\;
\frac{N}{K-N} = A\,e^{rt}, \quad A = \frac{N_0}{K - N_0}
$$

Λύνοντας ως προς $N$ καταλήγουμε στον **κλειστό τύπο της λογιστικής λύσης**:

$$
\boxed{\;N(t) \;=\; \dfrac{K}{1 + \left(\dfrac{K - N_0}{N_0}\right) e^{-rt}}\;}
$$

**Έλεγχος** (sanity checks):

- Στο $t = 0$: $N(0) = K \big/ \left(1 + \tfrac{K-N_0}{N_0}\right) = K \cdot \tfrac{N_0}{K} = N_0$. ✓
- Στο $t \to \infty$: το $e^{-rt} \to 0$, άρα $N \to K$. ✓
- Αν $N_0 \ll K$, για μικρά $t$ ο όρος $\tfrac{K-N_0}{N_0} \approx \tfrac{K}{N_0}$ είναι τεράστιος, και η συμπεριφορά είναι σχεδόν εκθετική (γιατί;). ✓

---

# R: αριθμητική λύση με `deSolve`

Τώρα που ξέρουμε τι περιμένουμε, ας προσομοιώσουμε:

```{r logistic-ode}
logistic_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N * (1 - N / K)
    list(dN)
  })
}

parameters <- c(r = 0.1, K = 10000)
state      <- c(N = 100)
times      <- seq(0, 200, by = 0.5)

out_log <- as.data.frame(
  ode(y = state, times = times, func = logistic_model, parms = parameters)
)

head(out_log)
```

## Σύγκριση αναλυτικής και αριθμητικής λύσης

```{r compare-analytical-numerical}
N0 <- 100
K  <- 10000
r  <- 0.1

logistic_analytical <- function(t, N0, r, K) {
  K / (1 + ((K - N0) / N0) * exp(-r * t))
}

df_log <- data.frame(
  t          = times,
  Αναλυτική  = logistic_analytical(times, N0, r, K),
  Αριθμητική = out_log$N
)

ggplot(df_log, aes(x = t)) +
  geom_line(aes(y = Αναλυτική),  linewidth = 1.2, colour = "steelblue") +
  geom_point(aes(y = Αριθμητική),
             data = df_log[seq(1, nrow(df_log), by = 25), ],
             colour = "tomato", size = 1.8) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95, y = K * 1.04,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Λογιστική αύξηση: αναλυτική (γραμμή) vs αριθμητική (κουκίδες)",
    x = "Χρόνος (ημέρες)",
    y = "Πληθυσμός N(t)"
  )

cat("Μέγιστο σφάλμα:", format(max(abs(df_log$Αναλυτική - df_log$Αριθμητική)),
                              scientific = TRUE), "\n")
```

---

# Η σιγμοειδής καμπύλη και τα τρία στάδια

Παρατηρήστε προσεκτικά τη χρονική τροχιά. Έχει τρία διακριτά στάδια:

```{r three-stages}
stages <- data.frame(
  t = times,
  N = logistic_analytical(times, N0, r, K)
)

# Σημείο καμπής: N = K/2
t_inflection <- (1/r) * log((K - N0) / N0)

ggplot(stages, aes(x = t, y = N)) +
  geom_line(linewidth = 1.1, colour = "steelblue") +
  geom_hline(yintercept = K,   linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = K/2, linetype = "dotted", colour = "grey50") +
  geom_vline(xintercept = t_inflection, linetype = "dotted", colour = "tomato") +
  annotate("point", x = t_inflection, y = K/2, size = 3, colour = "tomato") +
  annotate("text",  x = t_inflection + 5, y = K/2,
           label = "Σημείο καμπής (N=K/2)\nμέγιστος ρυθμός", hjust = 0) +
  annotate("text", x =  20, y = 1500, label = "1. Σχεδόν εκθετικό",     hjust = 0) +
  annotate("text", x =  60, y = 5500, label = "2. Μεταβατική φάση\n(επιτάχυνση → επιβράδυνση)", hjust = 0) +
  annotate("text", x = 130, y = 9300, label = "3. Κορεσμός προς K",       hjust = 0) +
  labs(
    title = "Η σιγμοειδής (S-shaped) λογιστική καμπύλη",
    x = "Χρόνος (ημέρες)",
    y = "N(t)"
  )

cat("Σημείο καμπής στο t =", round(t_inflection, 1), "ημέρες (όπου N = K/2 =", K/2, ")\n")
```

1. **Σχεδόν εκθετική φάση** (μικρά $N$): ο όρος $1 - N/K$ είναι κοντά στο 1, άρα $\frac{dN}{dt} \approx rN$.
2. **Σημείο καμπής** ($N = K/2$): εδώ ο ρυθμός $\frac{dN}{dt}$ μεγιστοποιείται.
3. **Φάση κορεσμού** ($N \to K$): ο όρος $1 - N/K \to 0$, άρα ο ρυθμός φθίνει εκθετικά προς το μηδέν.

---

# Σύγκριση: εκθετικό vs λογιστικό

Με ίδιο $r$, ίδιο $N_0$, αλλά διαφορετικό $K$:

```{r exp-vs-log}
df_compare <- data.frame(
  t          = times,
  Εκθετικό   = N0 * exp(r * times),
  Λογιστικό  = logistic_analytical(times, N0, r, K)
)

# Long format για ggplot
df_long <- data.frame(
  t       = rep(times, 2),
  N       = c(df_compare$Εκθετικό, df_compare$Λογιστικό),
  Μοντέλο = factor(rep(c("Εκθετικό", "Λογιστικό"), each = length(times)))
)

ggplot(df_long, aes(x = t, y = N, colour = Μοντέλο)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  scale_y_log10() +
  annotate("text", x = max(times) * 0.95, y = K * 1.3,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Εκθετική vs λογιστική αύξηση (λογαριθμικός άξονας)",
    subtitle = "Στις πρώτες ημέρες είναι σχεδόν ίδιες — μετά αποκλίνουν δραματικά",
    x = "Χρόνος (ημέρες)",
    y = "N(t) [λογαριθμική κλίμακα]"
  ) +
  theme(legend.position = "bottom")
```

Στις πρώτες $\approx 30$ ημέρες οι δύο καμπύλες είναι **σχεδόν αξεχώριστες**. Αυτό μας δίνει μια προειδοποίηση: στα πρώιμα στάδια ενός πληθυσμιακού φαινομένου, **δεν μπορούμε εύκολα να ξεχωρίσουμε** αν αυξάνεται εκθετικά ή λογιστικά — και τα δύο μοντέλα ταιριάζουν εξίσου καλά. Η διαφορά αναδύεται μόνο όταν πλησιάζουμε στο $K$.

---

# Τι συμβαίνει αν ξεκινήσουμε από **πάνω** από το $K$;

```{r overshoot}
times2 <- seq(0, 100, by = 0.5)

scenarios <- data.frame(
  t   = rep(times2, 3),
  N   = c(logistic_analytical(times2, N0 = 100,   r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 5000,  r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 15000, r = 0.1, K = 10000)),
  N0  = factor(rep(c("N0 = 100",  "N0 = 5.000",  "N0 = 15.000 (overshoot)"),
                   each = length(times2)))
)

ggplot(scenarios, aes(x = t, y = N, colour = N0)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  annotate("text", x = max(times2) * 0.95, y = K * 1.04,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Λογιστική αύξηση: τρία διαφορετικά αρχικά σενάρια",
    subtitle = "Όλα συγκλίνουν στο K, ανεξάρτητα από την αρχική συνθήκη",
    x = "Χρόνος (ημέρες)",
    y = "N(t)"
  ) +
  theme(legend.position = "bottom")
```

Παρατηρήστε ότι η περίπτωση $N_0 = 15.000 > K$ **μειώνεται** προς το $K$. Αυτό βγαίνει αυτόματα από την εξίσωση: αν $N > K$, τότε $1 - N/K < 0$, άρα $\frac{dN}{dt} < 0$.

Αυτό μοντελοποιεί π.χ. την περίπτωση που εισάγουμε υπερβολικό αριθμό νούφαρων σε μια λίμνη — ο πληθυσμός θα **πέσει** στο φέρον επίπεδο.

---

# Εκτίμηση παραμέτρων από δεδομένα

Στην πράξη οι βιολόγοι έχουν δεδομένα και θέλουν να **εκτιμήσουν** τα $r$ και $K$. Αυτό γίνεται με μη-γραμμική παλινδρόμηση (`nls()` στην R).

```{r fit-logistic}
# Προσομοιώνουμε δεδομένα με θόρυβο
set.seed(123)
sample_days <- seq(0, 150, by = 10)
true_N      <- logistic_analytical(sample_days, N0 = 100, r = 0.1, K = 10000)
observed_N  <- true_N * exp(rnorm(length(sample_days), 0, 0.08))

# Προσαρμογή του λογιστικού μοντέλου
fit <- nls(observed_N ~ K / (1 + ((K - N0)/N0) * exp(-r * sample_days)),
           start = list(K = 8000, r = 0.05, N0 = 200))

summary(fit)$coefficients
```

```{r plot-fit}
fit_curve <- data.frame(
  t = times,
  N = predict(fit, newdata = data.frame(sample_days = times))
)

ggplot() +
  geom_point(data = data.frame(t = sample_days, N = observed_N),
             aes(x = t, y = N), size = 2.5, colour = "tomato") +
  geom_line(data = fit_curve, aes(x = t, y = N),
            linewidth = 1, colour = "steelblue") +
  geom_hline(yintercept = coef(fit)["K"],
             linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95,
           y = coef(fit)["K"] * 1.05,
           label = paste0("Εκτίμηση K = ", round(coef(fit)["K"])),
           colour = "grey30", hjust = 1) +
  labs(
    title = "Προσαρμογή λογιστικού μοντέλου σε δεδομένα",
    subtitle = paste0("Εκτιμήσεις: r = ", round(coef(fit)["r"], 3),
                      ", K = ",  round(coef(fit)["K"]),
                      ", N0 = ", round(coef(fit)["N0"])),
    x = "Χρόνος (ημέρες)",
    y = "Πληθυσμός N(t)"
  )
```

> **Πρακτική σημασία:** Τα `nls()` εκτιμώμενα $\hat{r}$ και $\hat{K}$ έρχονται με τυπικό σφάλμα — άρα έχουμε στατιστικά διαστήματα εμπιστοσύνης για τις παραμέτρους. Αυτό είναι κρίσιμο για να συγκρίνουμε δύο πληθυσμούς ή δύο πειραματικές συνθήκες.

---

# Πέρα από τη λογιστική: άλλες λειτουργικές μορφές

Όπως αναφέραμε, η γραμμική μορφή $r_{\text{eff}}(N) = r(1 - N/K)$ είναι **μία επιλογή**. Δύο σημαντικές εναλλακτικές:

## θ-λογιστική

$$
\frac{dN}{dt} = r N \left(1 - \left(\frac{N}{K}\right)^{\theta}\right)
$$

- $\theta = 1$ → κλασική λογιστική.
- $\theta > 1$ → ο πληθυσμός μένει κοντά στο $K$ για περισσότερο χρόνο πριν επιβραδύνει (π.χ. $r$-στρατηγικές).
- $\theta < 1$ → η επιβράδυνση ξεκινά νωρίς (π.χ. ευαίσθητοι σε ανταγωνισμό).

## Gompertz

$$
\frac{dN}{dt} = r N \ln\!\left(\frac{K}{N}\right)
$$

Χρησιμοποιείται κυρίως σε αύξηση όγκων (π.χ. καρκινικοί όγκοι, ανάπτυξη οργανισμών). Επιβραδύνει **νωρίτερα** από τη λογιστική.

```{r alt-models}
gompertz <- function(t, N0, r, K) {
  K * exp(log(N0/K) * exp(-r * t))
}
theta_logistic <- function(t, N0, r, K, theta) {
  # Αριθμητική λύση
  parms <- c(r = r, K = K, theta = theta)
  func  <- function(t, state, parms) {
    with(as.list(c(state, parms)),
         list(r * N * (1 - (N/K)^theta)))
  }
  out <- ode(y = c(N = N0), times = t, func = func, parms = parms)
  as.data.frame(out)$N
}

t_vec <- seq(0, 200, by = 0.5)
df_alt <- data.frame(
  t        = rep(t_vec, 4),
  N        = c(logistic_analytical(t_vec, 100, 0.1, 10000),
               gompertz(t_vec,             100, 0.1, 10000),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 0.5),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 3)),
  Μοντέλο  = factor(rep(c("Λογιστικό (θ=1)",
                          "Gompertz",
                          "θ-λογιστικό (θ=0.5)",
                          "θ-λογιστικό (θ=3)"),
                        each = length(t_vec)),
                    levels = c("Λογιστικό (θ=1)",
                               "θ-λογιστικό (θ=0.5)",
                               "θ-λογιστικό (θ=3)",
                               "Gompertz"))
)

ggplot(df_alt, aes(x = t, y = N, colour = Μοντέλο)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 10000, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Διαφορετικές μοντελοποιήσεις του κορεσμού",
    subtitle = "Όλα τα μοντέλα έχουν ίδιο r=0.1, K=10.000, N0=100",
    x = "Χρόνος (ημέρες)", y = "N(t)"
  ) +
  theme(legend.position = "bottom")
```

> **Παιδαγωγικό σχόλιο:** Το ότι όλα αυτά τα μοντέλα δίνουν *διαφορετικές* καμπύλες είναι το πιο σημαντικό μήνυμα. Σε πραγματικά δεδομένα η επιλογή του μοντέλου **δεν προκύπτει αυτόματα** — εξαρτάται από τη βιολογία του οργανισμού και πρέπει να ελεγχθεί στατιστικά (π.χ. AIC).

---

# Ασκήσεις

1. **Ευαισθησία στο $K$.** Σχεδιάστε λογιστικές καμπύλες με $r=0{,}1$, $N_0=100$ και $K \in \{1.000,\ 5.000,\ 10.000,\ 50.000\}$ στο ίδιο διάγραμμα. Τι αλλάζει στη μορφή της καμπύλης;

2. **Χρόνος για να φτάσουμε στο $90\%$ του $K$.** Λύστε αναλυτικά: για ποιο $t$ ισχύει $N(t) = 0{,}9 K$; Επαληθεύστε αριθμητικά.

3. **Σημείο καμπής.** Δείξτε ότι ο ρυθμός $\frac{dN}{dt}$ μεγιστοποιείται στο $N = K/2$. (*Υπόδειξη:* παραγωγίστε τον $\frac{dN}{dt}$ ως προς $N$ και βρείτε τη ρίζα.)

4. **Από τα δεδομένα στο $r$.** Σας δίνονται μετρήσεις:
   ```r
   day <- c(0, 10, 20, 40, 60, 90, 120, 180)
   N   <- c(120, 280, 620, 2400, 5800, 9100, 9750, 9970)
   ```
   Προσαρμόστε το λογιστικό μοντέλο με `nls()` και υπολογίστε τα διαστήματα εμπιστοσύνης για τα $\hat{r}$ και $\hat{K}$ (`confint(fit)`).

5. **Εννοιολογική.** Στη συζήτηση αναφέραμε *διάφορους* μηχανισμούς που μπορούν να δικαιολογήσουν τη μείωση του κατά κεφαλήν ρυθμού (χώρος, φως, θρεπτικά, παθογόνα). Για τα νούφαρα στη λίμνη μας, ποιον μηχανισμό θεωρείτε **πιο σημαντικό** και γιατί; Πώς θα τον επαληθεύατε πειραματικά;

6. **Πέρα από τη μέση τιμή.** Όλα τα μοντέλα εδώ είναι ντετερμινιστικά — δεν υπάρχει τύχη. Σε πραγματικούς πληθυσμούς όμως υπάρχει στοχαστικότητα (γεννήσεις/θάνατοι είναι τυχαία γεγονότα). Δοκιμάστε να προσομοιώσετε στοχαστική λογιστική αύξηση: σε κάθε χρονικό βήμα, προσθέστε στον $\frac{dN}{dt}$ έναν τυχαίο όρο `rnorm(1, 0, σ)`. Πώς αλλάζει η συμπεριφορά για διαφορετικά $\sigma$;

---

# Επόμενο βήμα

Μέχρι τώρα έχουμε εξετάσει έναν *μοναδικό* πληθυσμό. Στην πραγματικότητα όμως οι πληθυσμοί **αλληλεπιδρούν** μεταξύ τους:

- **Θηρευτής–θήραμα**: εξισώσεις Lotka–Volterra
- **Ανταγωνισμός μεταξύ ειδών**: επεκτεταμένη λογιστική με δύο ή περισσότερες εξισώσεις
- **Συμβίωση**: αμοιβαία ωφέλιμες αλληλεπιδράσεις

Αυτά είναι **συστήματα διαφορικών εξισώσεων** και θα είναι το θέμα του επόμενου εργαστηρίου.

---

# Στοιχεία συνεδρίας R

```{r session-info}
sessionInfo()
```
